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Abstract 


Probability  of  close  approach  is  the  probability  that  two  satellites  will 
be  within  some  specified  distance  threshold  of  each  other  at  a  random  time 
within  a  specified  time  interval.  In  this  paper,  methods  were  developed  to 
calculate  probability  of  close  approach  between  two  satellites.  To  simplify 
the  analysis,  the  investigation  was  restricted  to  satellite  orbits  and  time 
intervals  where  the  mean  anomaly  of  both  satellites  can  be  treated  as 
independent,  uniformly  distributed  random  variables.  In  addition,  all  orbital 
parameters,  except  for  mean  anomaly,  were  assumed  to  be  constant  over 
time.  This  means  that  all  the  methods  developed  in  this  paper  to  calculate 
the  probability  of  close  approach  will  only  be  valid  over  very  long  time 
intervals  where  the  ratio  of  the  orbital  periods  of  the  two  satellites  can  be 
approximated  as  an  irrational  number.  Likewise,  there  can  be  no 
perturbations  in  the  orbital  parameters  of  both  satellites. 

The  first  method  developed  was  a  general  method  for  calculating  the 
probability  of  close  approach  between  two  satellites  in  elliptical  orbits.  The 
method  requires  numerical  integration  and  direct  solution  of  the  roots  of  a 
4th  order  polynomial  during  each  numerical  integration  step. 

Another  method  was  developed  for  calculating  the  probability  of  close 
approach  between  two  satellites  in  circular  orbits.  This  method  still  requires 
numerical  integration  to  obtain  a  solution,  but  in  this  case  a  direct  solution 
was  found  for  the  limits  of  integration.  Futhermore.  the  calculations 
required  during  each  numerical  integration  step  are  much  simpler  than  those 
required  to  calculate  the  probability  of  close  approach  with  elliptical  orbits. 

vii 


Finally,  a  direct  solution  for  approximate  probability  of  close  approach 
between  two  satellites  in  circular  orbits  was  developed  for  the  case  where  the 
angle  between  the  orbital  planes  of  both  satellites  is  not  small  and  the 
probability  of  close  approach  is  small. 

Both  the  elliptical  orbit  and  the  circular  orbit  methods  of  computing 
probability  of  close  approach  yielded  results  that  compare  favorably  with 
estimates  of  probability  of  close  approach  derived  from  statistical  simulations. 


I.  Introduction 


There  are  a  variety  of  problems  where  the  close  approach  of  two 
satellites  is  of  interest.  Here,  close  approach  of  two  satellites  is  defined  as 
ocurring  whenever  the  distance  between  two  satellites  is  less  than  or  equal  to 
some  distance  threshold  dra.  When  the  position  and  velocity  of  both 

satellites  are  well  known  the  actual  time  and  duration  of  each  close  approach 
can  be  predicted.  However,  if  the  time  of  interest  cannot  be  predicted,  then 
a  deterministic  approach  to  the  close  approach  of  two  satellites  can  no  longer 
be  used. 

The  purpose  of  this  paper  is  to  develop  methods  to  calculate  the 
probability  of  close  approach  between  two  satellites  at  a  uniformly 
distributed  random  time  within  a  specified  time  interval.  To  simplify  the 
analysis,  the  investigation  is  restricted  to  satellite  orbits  and  time  intervals 
where  the  mean  anomaly  of  both  satellites  can  be  treated  as  independent, 
uniformly  distributed  random  variables.  In  addition,  all  other  orbital 
parameters  are  assumed  to  be  constant  over  time.  Because  of  these 
restrictions,  the  methods  developed  to  calculate  probability  of  clc'e  approach 
are  only  valid  over  very  long  time  intervals  with  some  restrictions  to  the 
ratio  of  the  orbital  periods  of  the  two  satellites  (see  the  Theory  section  of 
chapter  III). 

In  general,  the  goal  is  to  come  up  with  a  way  to  calculate  the 
probability  of  close  approach  between  two  satellites  in  elliptical  orbits.  This 
general  method  can  also  be  used  to  calculate  the  probabilit)  of  close 
approach  between  two  satellites  in  circular  orbits,  but  it  is  computationally 


1 


cheaper  to  use  a  method  designed  specifically  to  calculate  th*  probability  of 
close  approach  between  two  satellites  in  circular  orbits.  Similarly,  there  are 
some  special  cases  where  an  approximate  method  for  calculating  the 
probability  of  close  approach  gives  adequate  accuracy  at  much  less 
computational  expense.  For  these  reasons,  three  different  methods  will  be 
developed  to  calculate  probability  of  close  approach.  The  first  method  is  a 
general  method  for  calculating  the  probability  of  close  approach  between  two 
satellites  in  elliptical  orbits  (see  chapter  III).  The  second  method  is  for 
calculating  the  probability  of  close  approach  between  two  satellites  in  circular 
orbits  (see  chapter  IV).  The  last  method  is  for  calculating  the  approximate 
probability  of  close  approach  between  two  satellites  in  circular  orbits  where 
the  probability  of  close  approach  is  small  and  the  angle  between  the  orbital 
planes  of  the  two  satellites  is  not  small  (see  chapter  IV). 

Finally,  to  verify  that  the  three  methods  are  correct,  the  probability  of 
close  approach  will  be  computed  (using  the  three  methods,  where  applicable) 
for  a  variety  of  orbital  test  cases,  and  the  results  will  be  compared  to  values 
derived  from  statistical  simulations. 


II.  Background 


work  has  been  done  in  investigating  the  probability  of  collision 
between  orbiting  bodies  (3,  6,  and  8).  Probability  of  collision  is  typically 
defined  as  the  probability  that  one  orbiting  body /satellite  will  come  within 
some  distance  threshold  of  another  satellite  one  or  more  times  within  a 
specified  time  interval.  When  dealing  with  probability  of  satellite  collision, 
the  distance  threshold  used  is  typically  very  small,  since  it  is  directly  related 
to  the  physical  size  of  the  satellites  involved. 

More  recently,  work  has  been  done  in  investigating  the  probability  of 
satellite  intercept  between  satellites  in  circular  orbits  (9).  Here,  a  satellite 
intercept  is  defined  as  occurring  when  two  satellites  are  within  some  distance 
threshold  of  each  other  for  at  least  some  specified  length  of  time. 
Probability  of  satellite  intercept  is  the  probability  that  one  or  more  satellite 
intercepts  will  occur  within  some  specified  time  interval  that  begins  at  some 
uniformly  distributed  random  time  within  another  much  larger  time  interval. 

Probability  of  collision  and  probability  of  intercept  are  two  examples  of 
probabilistic  measures  dealing  with  satellite  proximity.  This  paper  introduces 
a  new  probabilistic  measure  of  satellite  proximity,  called  probability  of  close 
approach.  Probability  of  close  approach  is  the  expected  fraction  of  a 
specified  time  interval  over  which  the  distance  between  the  two  satellites  is 
less  than  or  equal  to  some  distance  threshold.  For  very  long  time  intervals, 
the  probability  of  close  approach  equals  the  sum  of  the  durations  of  all  the 
close  approaches  that  occur  within  the  specified  time  interval  divided  by  the 
length  of  the  specified  time  interval. 

Probability  of  close  approach  is  very  different  from  probability  of 


3 


collision.  When  the  probability  of  one  or  more  collisions  between  two 
satellites  within  a  very  long  time  interval  approaches  1.0,  the  computed 
probability  of  close  approach  can  approach  zero.  The  reason  for  this  is  that, 
regardless  of  the  number  of  collisions,  the  actual  time  that  two  satellites 
spend  within  the  collision  distance  threshold  can  be  a  very  small  fraction  of 
the  length  of  the  time  interval  of  interest. 

Probability  of  close  approach  is  closer  in  concept  to  probability  of 
intercept,  but  there  are  still  major  differences.  Probability  of  close  approach 
places  no  requirement  on  the  duration  of  the  close  approach,  and  close 
approaches  that  occur  after  the  uniformly  distributed  random  time  are  of  no 
interest.  Despite  these  differences,  probability  of  intercept  and  probability  of 
close  approach  share  three  major  assumptions.  First,  the  time  of  interest  is 
assumed  to  be  a  uniformally  distributed  random  time  within  some  very  long 
time  interval.  Second,  the  mean  anomalies  of  both  satellites  are  assumed  to 
be  independent,  uniformly  distributed  random  variables  (probability  of 
intercept  was  derived  only  for  circular  orbits,  where  mean  anomaly  always 
equals  true  anomaly).  Finally,  all  other  orbital  parameters  are  assumed  to 
be  constant  over  time. 

Probability  of  close  approach  is  different  from  probability  of  intercept, 
just  as  their  purposes  are  different.  When  it  is  important  that  one  or  more 
intercepts  occur  between  two  satellites,  all  within  a  specified  time  interval 
starting  at  some  random  time,  a  high  probability  of  intercept  is  desirable. 
When  it  is  important  that  one  satellite  spend  as  much  of  its  orbital  lifetime 
as  possible  within  some  arbitrary  distance  threshold  of  another  satellite,  then 
a  high  probability  of  close  approach  is  desirable. 
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III.  Probability  of  Close  Approach  Between 
Satellites  in  Elliptical  Orbits 

The  purpose  of  this  chapter  is  to  develop  a  method  for  calculating  the 
probability  of  close  approach  between  two  satellites  in  elliptical  orbits. 

Theory 

The  derivation  of  probability  of  close  approach  can  be  broken  up  into 
three  major  parts.  This  section  identifies  these  major  parts  and  describes 
how  they  can  be  put  together  to  calculate  probability  of  close  approach. 
The  next  three  sections  of  chapter  III  then  completes  the  solution  for  each 
of  the  three  parts. 

By  definition,  probability  of  close  approach  is  the  probability  that  the 
distance  between  two  satellites  will  be  less  than  some  distance  threshold  at  a 
uniformly  distributed  random  time  within  a  specified  time  interval.  To 
simplify  the  analysis,  two  basic  assumptions  were  made.  First,  all  orbital 
elements,  except  for  mean  anomaly,  are  assumed  to  be  constant  over  time. 
Second,  the  mean  anomalies  of  both  satellites  are  assumed  to  be  independent 
random  variables  that  are  uniformly  distributed  between  0  and  2n.  The 
first  assumption  is  valid  when  there  are  no  perturbations  to  the  orbital 
elements  of  the  two  satellites.  When  is  the  second  assumption  valid?  At  a 
uniformly  distributed  time  within  a  specified  time  interval,  the  mean  anomaly 
of  both  satellites  can  be  represented  by  (1:33,  185) 

M,  =  (  2tt  t,Pl  x  Mto  )  MOD  2 n  (1) 

M,  =  (  2rr  t/P,  4  )  MOD  2tt  (2) 
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where 


t  =  time  from  start  of  time  interval  of  interest. 

=  mean  anomaly  of  satellite  1. 

Mlo  =  mean  anomaly  of  satellite  1  at  t  =  0. 

=  mean  anomaly  of  satellite  2. 

Mto  =  mean  anomaly  of  satellite  2  at  t  =  0. 

Pj  =  orbital  period  of  satellite  1. 

P,  =  orbital  period  of  satellite  2. 

and  the  general  function  X  MOD  Y  represents  the  remainder  of  X  divided 
by  Y.  For  purposes  of  this  analysis,  t  is  a  random  variable  uniformly 
distributed  between  0  and  the  duration  of  the  time  interval  of  interest. 

Since  mean  anomaly  is  a  linear  function  of  t  (see  Eqs  (1)  and  (2)),  Mj  and 

M,  are  also  uniformly  distributed  random  variables  when  the  duration  of  the 

time  interval  is  less  than  the  orbital  period  of  both  satellites,  or  when  the 
duration  of  the  time  interval  is  equal  to  some  integer  multiple  of  the  period 
of  both  satellites.  Furthermore,  over  very  long  time  intervals  (over  100 
orbital  periods)  and  M,  approximate  (within  1%)  random  variables  that 

are  uniformly  distributed  between  0  and  2n.  Therefore,  over  long  time 
intervals  the  mean  anomalies  of  both  satellites  can  be  treated  as  uniformly 
distributed  random  variables  (not  necessarily  independent). 

When  can  Mt  and  M?  be  considered  independent?  Let  to  be  some 

arbitrary  time  within  the  time  interval  of  interest,  and  let  n  be  some 

nonnegative  integer.  Using  Eqs  (1)  and  (2),  at  t  «  to  *f  nP,  the  mean 

anomalies  of  both  satellites  can  be  represented  by  the  following: 
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(3) 


“i  =  Mlt 

M2  =  [  2rr  (Pj/PjJn  4  Mw  ]  MOD  2tt  (4) 

where 

Mu  =  the  mean  anomaly  of  satellite  1  at  t  =  to. 

M^  =  the  mean  anomaly  of  satellite  2  at  t  =  to. 

Similarly,  at  t  =  t  4-  nP2  the  mean  anomalies  of  both  satellites  can  be 

represented  by 

Mj  =  [  2tt  (Pj/PjJn  4  Mlt  ]  MOD  2t r  (5) 

M2  =  U7%  (6) 

Mj  and  M}  can  be  considered  independent  so  long  as  the  results  of  Eqs  (4) 

and  (5)  are  uniformly  distributed  between  0  and  2 n.  Once  again,  let  n  be 
some  nonnegative  integer.  Also  let  nmax  be  the  maximum  number  of  orbital 
periods  within  the  time  interval  of  interest.  There  are  two  cases  in  which 

Mj  and  M2  can  be  considered  independent.  The  first  case  is  where  P/Pj 
and  P2/Pt  are  irrational,  and  for  n  between  0  and  nmax,  the  distances  from 
(Pi/P2)n  and  (P2/Pj)n  to  the  nearest  integer  are  not  less  than  one  divided 
by  nmax.  For  practical  purposes,  the  ratios  Pj/P,  and  P,/Pj  can  be 
considered  irrational  when  (Pj/P2)n  and  (P2/Pj)n  do  not  equal  integers  for  n 
between  0  and  nmax.  The  second  case  is  where  (Pj/PjJn  and  (Pt  P2)n 
equal  integers  for  some  value  of  n  less  than  nmax,  where  n  is  large  (lOOU-r). 
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and  either  nmax  MOD  n  is  large  (1000+),  or  the  distance  from  nmax 
divided  by  n  to  the  nearest  integer  is  small.  For  members  of  the  second 
case  where  nmax  MOD  n  is  large,  the  orbital  ratios  must  also  meet  the 
criteria  of  the  first  case  for  n  between  0  and  nmax  MOD  n. 

Generally  speaking,  if  Pj/P,  and  Pj/P,  a.*  irrational,  and  the  time 

interval  of  interest  is  very  long  (1000-1-  orbital  periods),  then  the  mean 
anomalies  of  both  satellites  can  be  considered  independent,  uniformly 
distributed  random  variables. 

As  discussed  above,  the  mean  anomalies  of  both  satellites  are  assumed 
to  be  independent  random  variables  that  are  uniformly  distributed  between 
of  0  to  2 it.  This  means  that  the  probability  density  functions  of  the  mean 
anomalies  of  the  two  satellites  can  be  represented  by  (5:72-73) 


AM,)  = 

-  1 
•  2t r 

for  0  £  M, 

<  2rr 

(7) 

AM,)  = 

l 

'  2 W 

for  0  £  M, 

<  2  n 

(8) 

where 

Mj  =  the  mean  anomaly  of  satellite  1. 

M3  =  the  mean  anomaly  of  satellite  2. 

Likewise,  the  joint  probability  density  function  of  the  mean  anomalies  of 
both  satellites  can  be  represented  by  (5:135,  139-140) 

f[ M..MJ  =  -I,  for  0  S  M  <  2 it  (9) 

'  1  *  4rr-  1 

0SM,<  2rr 
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The  probability  of  close  approach  (PCA)  can  be  computed  by 
integrating  the  joint  density  function  over  the  region  of  and  where 

the  distance  between  the  satellites  (d)  is  less  than  or  equal  to  some  distance 
threshold  dra.  Eq  (10)  is  the  formula  for  probability  of  close  approach 

between  the  two  satellites: 


S, 


(10) 


where  Rj  is  the  region  of  M,  over  which  a  close  approach  occurs,  given  Mr 
and  Rt  is  the  region  of  over  which  some  close  approach  with  satellite  2 
is  possible. 

To  simplify  analysis,  four  functions  will  be  defined.  At  this  point, 
these  functions  are  strictly  symbolic,  and  no  solution  for  these  functions 
exist.  The  four  functions  are  MjO/j),  Mjft',),  d(MltMt),  and  AMa(Mj). 

M^i/j)  is  the  mean  anomaly  of  satellite  1  as  a  function  of  the  true  anomaly 
of  satellite  1.  Similarly,  M,(i/t)  is  the  mean  anomaly  of  satellite  2  as  a 
function  of  the  true  anomaly  of  satellite  2.  The  function  d(MpMt) 
represents  the  distance  between  the  two  satellites  as  a  function  of  their  mean 
anomalies.  The  last  function  is  AMt(M(),  which  is  defined  as 


AM,(Mt) 
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Substituting  Eq  (11)  into  Eq  (10)  results  in 


PCA  =  f  (12) 

J  4jt* 

R, 

Given  Mr  the  distance  between  the  two  satellites  is  purely  a  function  of  M3 
(all  other  orbital  elements  are  assumed  constant).  Each  solution  to  the 
equation  d(MrM3)  =  dafH  will  be  referred  to  as  a  mean  anomaly  close 

approach  boundary.  The  reason  for  this  is  that  they  delimit  the  mean 
anomaly  regions  of  close  approach  between  both  satellites.  For  each  two  M, 

solutions  to  d(M1YMt)  =  dfH  ,  there  is  a  mean  anomaly  close  approach 
region  such  that 

d(MlfMt)  S  dra  for  M m  £  M,  £  M1/tJ 

^i/n  * 

where  Mt/n  is  used  to  represent  the  beginning  of  the  ifA  close  approach 
region,  and  is  used  to  represent  the  end  of  the  i th  close  approach 

region.  If  there  is  more  than  one  close  approach  region,  then  the  regions  are 
numbered  so  that  the  mean  anomaly  of  the  beginning  of  the  i+1  close 
approach  region  is  greater  than  or  equal  to  the  mean  anomaly  of  the  end  of 
the  i th  close  approach  region.  The  possible  range  of  each  mean  anomaly 

10 


close  approach  boundary  is  from  0  to  2rr.  The  exception  to  this  is  when  a 
close  approach  region  crosses  2tt.  In  this  case,  Ms/U  (the  beginning  of  the 

first  close  approach  region)  can  range  from  -2rr  to  0  so  that  the  close 
approach  region  that  crosses  27T  does  not  have  to  be  broken  into  two  parts. 
Given  all  this,  when  there  are  n  close  approach  regions  and  n  is  greater 
than  zero,  the  general  solution  for  is 

=  t  <  Mj/b  -  MJ/U  )  (12) 

tal 

where 

-2rr  £  MJ/U  £  2jt 

0  £  M,/(l  £  2rr  for  i  >  1 
0  £  M,/u  £  2rr 

d(M,,M,)  £  dra  for  M1/n  £  M,  £  MJ/tt 

d(M.,K../U)  = 
d(M,.Mw)  = 

^t/U  ^  ^l/tt 

MJ/J}  <  Mj/3,  for  n  >  1  and  j  >  i 

For  example,  when  there  are  four  Ma  solutions  to  d(Mt,Mj)  =  d^  ,  the 
solution  to  AM2(Mj)  can  be  represented  by  (also  see  Figure  1) 

AMj(Mj)  =  MJ/}J  -  MJ/ai  4  Mt/|J  -  MJ/n  (13) 
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where 


i 

B 

I 

I 

K 

I 

8 


-2n  £  M}/11  £  2n 

0  £  Mj/h  £  2n 

0  £  Mt/tl  £  2rr 

0  £  M}/11  £  2n 

d(MltMJ/u)  =  dra 

d(M,,MJ/11)  =  dra 

=  dffi 

d(M1,M,/n)  =  d„ 

^1/U  ^  ^I/il 

^  ^*/Jl 

^»/»i  ^  ^j/ii 

d(M,,M1)  £  d^ 

for 

^i/n  &  £  ^j/n 

d(M,,M,)  £  dTH 

for 

Mj/ji  *  M,  £  M„„ 

When  there  are  no  M,  solutions  to  d(Mt,Mt)  =  dra  ,  AM^Mj)  can  possess 
one  of  two  possible  values.  If  dfM^M,)  >  dm  for  0  £  M,  £  2tt  ,  then 
AMt(Mt)  must  equal  sero.  If  d(Mt>Mt)  £  dra  for  0  £  M,  £  2n  ,  then 
AM,(M,)  must  equal  2rr.  However,  even  when  there  are  no  M,  solutions  to 
d(M|tMt)  =  d^  ,  Eq  (12)  can  still  be  used  to  calculate  AM^Mj)  through 
the  following  procedure: 

If  dfM^M,)  >  dTO  for  0  £  M,  £  2n  (14) 

n  =  1 


u\sT*.  t’G'i  rvri.  fttWCiT  xFJfF, 


If  d(UvMt)  *  <*th  for  0  £  M,  £  2n 


(15) 


n  =  1 


Figure  1.  Description  of  Close  Approach  Boundaries 


With  non-circular  orbits,  true  anomaly  is  much  easier  to  work  with 
than  mean  anomaly.  With  this  in  mind,  let  d  now  be  represented  as  a 
function  of  true  anomaly  instead  of  mean  anomaly.  In  other  words,  d(i/t,t/2) 

now  represents  the  distance  between  the  two  satellites  as  a  function  of  the 
true  anomalies  of  both  satellites  (i/j  and  u%).  Given  the  true  anomaly  of 


,  V*  ^  *  v%  JV  .’i  urn 


satellite  1  the  distance  between  the  two  satellites  is  purely  a  function 

of  the  true  anomaly  of  satellite  2  (i/,).  Each  vt  solution  to  the  equation 
d(vvv2)  =  dra  will  be  referred  to  as  true  anomaly  close  approach  boundary. 
For  each  two  i/,  solutions  to  d(i/lti/,)  =  dra  ,  there  is  a  true  anomaly  close 
approach  region  such  that 


d(i/i,i/,)  £  dra  for  v 


'l/tt 


vm  *  vm 


where  i/J/u  is  used  to  represent  the  true  anomaly  of  the  beginning  of  the  itA 
close  approach  region,  and  represents  the  true  anomaly  of  the  end  of 

the  itA  close  approach  region.  Mean  anomaly  close  approach  boundaries  and 
true  anomaly  close  approach  boundaries  are  related  in  the  following  way: 


^i/li 

(16) 

=  M,(i/J/U) 

(17) 

where  M,(i/,)  is  the  mean  anomaly  of  satellite  2  as  a  function  of  true 
anomaly. 

Given  vv  AM,  can  now  be  expressed  as  a  function  of  the  true  anomaly 

of  satellite  1.  When  there  are  n  close  approach  regions  and  n  is  greater 
than  zero,  the  general  solution  for  AM,(i/|)  is  (also  see  Eq  (12)) 


=  t  (  M, (",,«)  -  M,(t/  )  ]  (18) 
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vLy-V  /  vv. 


where 


8 


a 


i 


s 


-2rr  £  $  2ir 

0  £  vw  £  2rr 

for  i  >  1 

0  S  ^  S  2rr 

*<V'.)  *  dTH 

for  i/|/u  i/1/B 

^(^It^j/u)  =  ^TH 

v\l\\  ^  vipa 

vm  *  "vn 

for  n  >  1  and  j  >  i 

When  there  are  no  solutions  to  i(vvv%)  =  ,  Bq  (IS)  can  still  he 

used  to  compute  AMt(i/t)  through  the  following  procedure  (also  see  Eqs  (14) 
and  (15)): 


If  d(t/|li's)  >  d^  for  0  »  i/f  £  2n 


(19) 


n  =  1 


t/*/n  ■ 0 


^l/M  ~  ® 


If  dU'j.i/j)  £  d^  for  0  £  v%  £  2n 


(20) 


n  -  I 


Vvn  =  u 
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By  changing  the  integration  variable  in  Gq  (12)  from  mean  anomaly  of 
satellite  1  to  true  anomaly  of  satellite  1T  the  equation  for  probability  of  close 
approach  becomes  (4:212) 

PCA  =  AMjfi/j) 

T. 

where  T}  is  the  region  of  vx  over  which  a  close  approach  is  possible. 

Using  Eq  (21)  to  compute  PCA  requires  integration  over  the  region  of 
vx  where  some  close  approach  with  satellite  2  is  possible.  Since  AM](t/|)  can 
be  discontinuous  for  some  i/)t  this  requires  that  the  limits  of  integration  be 
found  over  which  the  function  AMJ(i/|)  is  continuous  before  Eq  (21)  can  be 

integrated  analytically.  Unfortunately,  there  is  normally  no  closed  form 
solution  for  the  limits  of  integration.  Therefore,  to  calculate  probability  of 
close  approach,  Eq  (21)  must  be  numerically  integrated  over  the  complete  2n 
range  of  vy 

The  final  equation  for  probability  of  close  approach  between  two 
satellites  in  elliptical  orbits  is 


^  dl/.  i 


in 


dl/. 


(21) 
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where  for  each  numerical  integration  step,  AM2(i/j)  is  computed  using  Eqs 
(18)  -  (20). 

Three  more  things  are  needed  to  complete  the  solution  for  probability 
of  close  approach: 

1.  An  equation  for  M2(l/2)  (see  Eq  (18)) 

2.  A.i  equation  for  dM^i/jJ/dt/j  (see  Eq  (22)) 

3.  Given  and  dTH,  a  method  to  determine  the  close  approach 
boundaries  (i/2^u  and  i/2/l2)  (see  Eq  (18)  -  (20)) 

The  solutions  to  these  problems  are  the  subject  of  the  next  three  sections  of 
chapter  HI. 

Mean  Anomaly  as  a  Function  of  True  Anomaly 

The  purpose  of  this  section  is  to  derive  an  equation  for  mean 
as  a  function  of  true  anomaly. 

Eqs  (23)  and  (24)  are  the  well  known  equations  relating  true 
to  eccentric  anomaly  (2:62),  and  mean  anomaly  to  eccentric  anomaly 

tan(i//2)  =  {(l+e)/(l-e)]1'2  tan(E/2) 

M  =  E  -  e  sinE 


anomaly 

anomaly 

(1:85): 

(23) 

(24) 


where 

v  -  true  anomaly 
e  -  eccentricity 
E  -  eccentric  anomaly 
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Substituting 


fi  =  [(l-e)/(l+e)]1'*  (25) 

into  Eq  (23),  and  then  rewriting  it  as  an  equation  for  eccentric  anomaly  as 
a  function  of  true  anomaly,  results  in 

E  =  2  tan”l[/?  tan(i//2)]  (26) 

Substituting  Eq  (26)  into  Eq  (24),  yields  an  equation  for  mean  anomaly  as  a 
function  of  true  anomaly: 

M  =  2  tan"1^  tan(i//2)]  -  e  sin{2  tan‘l[£  tan(i//2)]}  (27) 


Substituting 


f  =  fi  tan(i//2)  (28) 

into  Eq  (27),  results  in  an  equation  for  mean  anomaly  as  a  function  of  i/: 

M  =  2  tan”1^  -  e  sin[2  tan“s^]  (29) 

By  applying  the  trigonometric  relationship  (7:190)  sin29  =  2  sinG  cosG  , 
where  0  =  tan~1^  ,  Eq  (29)  becomes 

M  =  2  tan~1^’  -  e  [2  sin(tan~1^)  cosftan*"1^)]  (30) 


Substituting  (7:193) 


sin(tan“^)  =  i/(\  +^),/J 


(31) 


cos(tan“1^)  =  1/(1+^*) l^2 


(32) 


into  Eq  (30),  results  in 

M  =  2  tan"1^  -  2e  [t/(l+fi)l^[l/(l+1'2)l,r\ 

M  =  2  [tan"1^  -  e  i//(  1+^2)]  (33) 

Substituting  Eq  (28)  back  into  Eq  (33),  yields 

M  =  2  {tan"1^  tan(z//2)]  -  e  0  tan(i//2)/[l  +  P  tan2(i//2)]}  (34) 

As  a  last  step,  multiply  the  right  half  of  Eq  (34)  by  cos2(i//2)/cos2(i//2)  and 
simplify: 

M ' 2  [  ]  (36) 

Eq  (35)  is  the  equation  for  mean  anomaly  as  a  function  of  true  anomaly, 
which  is  one  of  the  things  needed  to  compute  the  probability  of  close 
approach.  However,  when  using  Eq  (35)  to  compute  mean  anomaly  in  a 
computer  program,  the  program  should  first  check  the  value  of  v.  If  v 
equals  tt,  then  the  program  should  directly  set  mean  anomaly  to  n  instead 
of  trying  to  calculate  mean  anomaly  using  Eq  (35),  because  tan(rr/2)  is 
infinite.  Likewise,  if  v  equals  -rr,  then  the  program  should  directly  set 
mean  anomaly  to  ~rr.  When  v  is  not  equal  to  ±7T,  then  Eq  (35)  can  safely 
be  used  to  compute  the  mean  anomaly. 
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Derivative  of  Mean  Anomaly  With  Respect  to  True  Anomaly 


i 


M(y)  represents  mean  anomaly  as  a  function  of  true  anomaly  (see  Eq 
(35)).  The  purpose  of  this  section  is  to  solve  for  the  derivative  of  M(i/) 
with  respect  to  true  anomaly. 

Differentiating  Eq  (33)  yields 

dM  =  2  [  1/(1  +f2)  -  e/U-hf'*)  -  2  ty  (36) 


S! 

x 


8 


i 


8 

!  9 


JV 


5 


t  * 
i  t 


A 

\  * 


I  V. 

! 4 


t  if 


i  E 
i * 


By  using  (l-f^2)3  as  a  common  denominator  in  all  three  terms  above,  Eq 
(36)  can  be  simplified  into 


du  " 2  [  ‘“(ufy1*1  ]  d* 


Differentiating  Eq  (28),  yields 

df  =  (/?/2)  sec,(i//2)  di/ 

Substituting  Eqs  (28),  and  (38)  into  Eq  (37),  results  in 


jM  - 2  [  'T.yjwl m) 


Substituting 


(l+e)0s  =  1-e 
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(37) 


(38) 


(39) 


. . . . . . —  . . . 


into  Eq  (39),  results  in 


dM 


=  2 


dM  =  2 


'i;  t  * 

('[T"+  u3?y  ]  vm  ~*W2>  if 


(40) 


Using  the  trigonometric  relationship 

l-ftan2(i//2)  =  sec2(i//2) 
Eq  (40)  can  be  simplified  to 


dU/du  =  2 


(l-e)secV/2)  ' 
[1  +  (P  tanJ(i//2)]J  _ 


(0/2)  secVtz/2) 


dM/di/  = 


(l-e)fl  sec4(t//2) 
0‘  tan}(i//2)]J 


[1  + 


(41) 


As  a  last  step,  multiply  the  right  side  of  Eq  (41)  by  cos*(v/2)/cos*{v/2 )  , 

and  then  simplify: 


dM/dl/  fcos5(i//2)  (+  ^sin»(t//2)]» 


(42) 


This  removes  any  potential  numerical  problems  at  v  —  in. 

Eq  (42)  is  the  equation  for  the  derivative  of  mean  anomaly  with 
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respect  to  true  anomaly,  which  is  the  second  of  three  things  needed  to 
compute  the  probability  of  close  approach. 

Finding  Close  Approach  Boundaries  Given 
the  True  Anomaly  of  Satellite  1 

Given  vv  when  there  are  2n  v2  solutions  to  d(vvvj  =  dTH  and  n  is 

positive,  there  are  n  true  anomaly  close  approach  regions.  The  close 
approach  boundaries  of  the  ith  close  approach  region  can  be  represented  by 

vm  and  vw  where 


-2-n  g  vyn  £  2 n 


0  S  i/J/u  £  2n  for  i  >  1 


0  £  vm  £  2rr 


d(VyV})  £  dTH  for  i')ln  £  i/,  £  vi/a 

A(vvvm)  =  dTH 


d^x -«'*«)  =  d 


TH 


t/2/U  S  ^2/12 


^2/12  **  V2i\l 


for  n  >  1  and  j  >  i 


The  close  approach  boundaries  for  each  close  approach  region  arc  required  to 
calculate  AM2(l>2)  (see  Eq  (18)),  and  a  method  to  determine  the  close 

approach  boundaries  of  each  close  approach  region  is  the  last  thing  required 
to  complete  the  solution  for  probability  of  close  approach  between  two 
satellites  in  elliptical  orbits.  The  purpose  of  this  section  is  to  develop  a 


OO 


4.4. 


method  to  calculate  vifn  and  i/J/l2  for  each  close  approach  region  given  the 
true  anomaly  of  satellite  1. 

When  numerically  integrating  Eq  (22),  only  vt  and  i/}  are  known  at 

the  beginning  of  each  numerical  integration  step.  To  compute  the  distance 
between  satellite  1  and  2,  the  position  vectors  of  both  satellites  must  be 
determined  within  a  common  cartesian  coordinate  frame.  For  convenience, 
the  perifocal  frame  of  satellite  2  was  selected. 

The  position  vector  of  satellite  2  in  the  perifocal  frame  of  satellite  2 
can  be  represented  by  (1:72) 

=  { xa  y*  za  ) 

where 

x2  =  r,  cos(i/a) 
y2  =  r2  sin(i/2) 
z2  =  0 

rj  —  Pj  /  (H“CjC°s(i/})] 
and 

e2  =  eccentricity  of  satellite  2 
p2  =  semi-latus  rectum  of  satellite  2 
=  true  anomaly  of  satellite  2 
r2  =  the  magnitude  of  r3/rJ 

Likewise,  the  position  vector  of  satellite  1  in  the  perifocal  frame  of  satellite  1 
can  be  represented  by  (1:72) 


(43) 

(44) 

(45) 

(46) 

(47) 
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ti,tx  =  <  x»  y.  2.  >  («) 

where 


x,  =  r,  cos(i/,) 

(49) 

y,  =  r,  sin(i/,) 

(50) 

zi  =  o 

(51) 

r,  =  P,  /  (l+ejcosti/,)] 

(52) 

and 

ej  =  eccentricity  of  satellite  1 
p,  =  semi-latus  rectum  of  satellite  1 
i/j  ==  true  anomaly  of  satellite  1 
=  the  magnitude  of 

Transforming  the  coordinate  frame  of  Tj  from  the  perifocal  frame  of  satellite 

1  to  the  perifocal  frame  of  satellite  2  can  be  performed  in  two  steps.  The 
first  step  is  to  transform  r,  from  the  perifocal  frame  of  satellite  1  to  the 

earth  centered  inertial  reference  frame.  The  last  step  is  to  transform  from 

the  inertial  frame  to  the  perifocal  frame  to  the  perifocal  frame  of  satellite  2. 

The  transformation  from  the  perifocal  frame  to  the  inertial  frame  can 
be  done  by  multiplying  the  position  vector  in  the  perifocal  frame  by  the 
following  transformation  matrix  (1:82-83): 


^»i 

K  R*  K 

^11  ^11 


(53) 
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where 


Rn  =  cos  0  cos  cj  -  sin  0  sin  cj  cos  i 

(54) 

Rl3  =  -cos  0  cos  cj  -  sin  f)  sin  o  cos  i 

(55) 

Rlf  =  sin  n  sin  i 

(56) 

Rai  =  sin  Q  cos  cj  -f  cos  0  sin  cj  cos  i 

(57) 

Rj2  =  -sin  fl  sin  cj  +  cos  0  cos  cj  cos  i 

(58) 

R}J  =  -cos  0  sin  i 

(59) 

R^  -  sin  cj  sin  i 

(60) 

RJ2  =  cos  cj  sin  i 

(61) 

Rjj  =  cos  i 

(62) 

and 

i  =  orbital  inclinition 
cj  =  argument  of  perigee 

0  =  longitude  of  the  ascending  node 

Regardless  of  the  perifocal  plane  that  the  position  vector 

transformed  into,  the  position  vector  of  satellite  1  in  the 

unchanged.  This  means  that 

of  satellite  1 

inertial  frame 

-i/i  = 

(63) 

-i/i  =  Ziw 

(64) 

Tj/^j 

(65) 
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where 


inclination  of  satellite  1 
inclination  of  satellite  2 


argument  of  perigee  of  satellite  1 

argument  of  perigee  of  satellite  2 

longitude  of  the  ascending  node  of  satellite  1 

longitude  of  the  ascending  node  of  satellite  2 

the  position  vector  of  satellite  1  in  the  inertial  frame 


the  position  vecto;  of  satellite  1  in  the  perifocal  plane  of 
satellite  1 


the  position  vector  of  satellite  1  in  the  perifocal  plane  of 
satellite  2 


the  transformation  matrix  to  transform  from  the  perifocal 
frame  of  satellite  1  to  the  inertial  frame 


the  transformation  matrix  to  transform  from  the  perifocal 
frame  of  satellite  2  to  the  inertial  frame 


Multiplying  both  sides  of  Eq  (65)  by  gr^i^cj^n,)  results  in  an  equation  for 


the  position  vector  of  satellite  1  in  the  perifocal  frame  of  satellite  2: 


fiftt  =  5(», r1/M 


The  transformation  matrix  R  is  orthogonal  (1:79-83),  so 


RT(i,,Oj,nj)  =  «r0,) 


Substituting  Eq  (67)  into  Eq  (66),  yields  the  final  equation  for  r.  : 


RT(i|,w,,nj)  R(i f, 


Eqs  (43)  through  (62),  and  Eq  (68)  make  it  possible  to  compute  the 
position  vectors  of  both  satellites  within  the  perifocal  frame  of  satellite  2. 
Given  that 


=  { x. 

y,  *,  } 

(69) 

-r5/,J  =  <  xs 

y,  o  } 

(70) 

the  distance  between  satellite  1  and  satellite  2  can  be  represented  by 

d  =  I  (x}  -  x,)»  +  (y,  -  y,)»  +  z\  r  (71) 

Simplifying  Eq  (70)  further,  yields 

d  =  (  x’-2x,x}+xj  +  y|-2y,y,+y{  +  i]  )l/1 
d  =  (  x|+yj+z{  +  x*+y*  -  2x,Xj  -  2yiy,  )*/»  (72) 

Substituting  rj  for  xj+yj+zj  and  rj  for  xj+yj,  Eq  (72)  becomes 

d  =  (  rj  +  rj  -  2XjXj  -  2yiy,  )*  (73) 


Substituting  Eqs  (44)  -  (47)  into  Eq  (73)  results  in 


d  = 


,  Pl _ toftcnW 

1  [l+e:cos(i/:)ll/3  l-fe5cos(j/?) 


2ylp,sin(tO  1*/* 

■  ■  V  1  (74) 

l-fe:cos(i/2) 


After  squaring  both  sides,  subtracting  rj  from  both  sides,  and  then 
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multiplying  both  sides  by  [l+eJcos(i/J)|,1  Eq  (74)  becomes 

(d’-rfHl+ejCosCi/,)]’  =  p’  -  2x)pJcos(t/})[l+e}cos(i/J)] 

-  2y jPjSin  (i/,) [1  +e,cos(i/j)]  (75) 

Simplifying  Eq  (75)  further 

(d3-r3)(l+2e2co$(i/2)+e3cos3(i/2)]  -  pj  +  2x1p2cos(t/2)[l+e2cos(i/2)] 

=  -2y1pJsin(i/2)[l-fe1cos(i/J)]  (76) 

Now  Eq  (76)  can  be  expressed  as  a  polynomial  of  cos(i/f): 

A  cos*(i/3)  +  B  cos(i/2)  +  C  =  -2y1p2sin(i/2)[H-e2cos(i/1)]  (77) 

where 


A  =  +  2x,Pj3 

(78) 

B  =  2[e,(d,-rj)  +  x,pj 

(79) 

C  =  d’-r^-pj 

(80) 

Squaring  both  sides  again,  Eq  (77)  becomes 

A3cos4(i/})  4  2ABcos,(i/2)  4  (2AC4B*)cos*(i/2)  4  2BCcos(i/2)  4  C3 

=  4y3p3sin3(i/2)[l-f2e2cos(i/2)-fejcos,(^J)]  (81) 

Using  the  trigonometric  relationship  1-cosVj  -  sinV2  in  Eq  (70)  yieldr 

A3cos4(t*2)  -f  2ABcos,(t/J)  4  (2AC4B3)cos3(i'2)  4  2BCcos(t'2)  4  C3 


=  4y*p*[l  -cos*(i/,)]  [1  +2e,cos(i/,)+e’cos,(t/,)] 

AW(i/,)  +  2ABcos‘(i/,)  +  (2AC+B*)cos‘(i/,)  +  2600)8(1/,)  +  C* 


=  4y’p*(l+2e,cos(«/,)+(«*-l)cos,(t/,)-2e,cos'(i/,)-e}cos4(t/,)]  (82) 

Now  Eq  (82)  can  be  expressed  as  a  single  4th  order  polynomial: 

PjCos^i/j)  +  Ptcos*(vt)  +  PJco5,(i/J)  +  P4cos(i/,)  -f  Pj  =  0  (83) 

where 

K  =  4yjp*  (84) 

P,  =  A1  +  Kc*  (86) 

P,  =  2AB  4-  2Ke,  (86) 

P,  =  2AC  +  B*  +  K(l-eJ)  (87) 

P4  =  2CB  -  2Ke,  (88) 

P4  -  Ca  -  K  (89) 


The  roots  of  a  4th  order  polynomial  can  be  solved  for  directly 
(7:103-106).  This  means  that,  given  i/)f  Eq  (83)  can  be  used  to  solve  for 

all  possible  values  of  cos(i/t)  where  the  distance  between  the  two  satellites  is 
equal  to  d^. 

A  special  case  occurs  when  the  absolute  value  of  yt  approaches  tero. 


When  yx  equals  sero,  Eq  (83)  reduces  to  (see  Eqs  (77)  to  (80)) 


A  cos,(i/J)  +  B  cos(i/,)  4  C  =  0 


(90) 


where 


A  =  e,[(d’-r*)  +  2x,pJ 

(91) 

B  =  2[e,(d*-r*)  +  x,pj 

(92) 

C  =  d*-rj-p| 

(93) 

This  also  means  that,  when  yt  equals  zero,  £q  (83)  is  equal  to  the  square  of 
a  2nd  order  polynomial.  Theoretically,  when  equals  zero,  the  roots  of  the 
square  of  £q  (90)  are  the  same  as  two  copies  of  the  roots  of  Eq  (90). 

However,  in  practical  applications,  this  is  not  the  case.  When  using  1EEF 
double  precision  arithmetic,  the  direct  solution  of  the  4th  order  roots  of  the 
square  of  Eq  (90)  can  result  in  a  pair  of  complex  conjugate  roots  for  each 

real  root  of  Eq  (90),  where  the  real  component  of  each  pair  of  complex 

conjugate  roots  would  equal  one  of  the  real  roots  of  Eq  (90),  and  the 

imaginary  component  would  be  some  small  value  on  the  order  of  10~*. 

The  addition  of  any  imaginary  number  to  an  otherwise  valid  solution 
for  cos(t/f),  makes  that  solution  unusable.  Because  of  this,  when  the 
absolute  value  of  yt  is  small,  two  copies  of  the  roots  of  Eq  (90)  should  be 

used,  instead  of  directly  solving  for  the  4th  order  roots  of  Eq  (83). 

Under  some  conditions,  it  is  possible  that  both  P(  and  P2  in  Eq  (83) 

will  equal  zero.  For  example,  when  satellite  2  is  in  a  circular  orbit  (e, 
equab  zero),  both  PI  and  P2  in  Eq  (83)  are  always  equal  to  zero.  This  case 
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can  be  handled  by  checking  the  values  of  Pt  and  Pf,  before  solving  for  the 
4th  order  roots  of  Eq  (83).  If  both  Px  and  P2  in  Eq  (83)  are  equal  to  zero, 

then  the  desired  solutions  are  the  roots  of  the  remaining  2nd  order  equation. 
Another  special  case  occurs  when  yx  approaches  zero  and  satellite  2  is  in  a 
circular  orbit.  Ignoring  the  fact  that  satellite  2  is  in  a  circular  orbit,  since 
y2  is  approximately  zero,  the  desired  solutions  should  be  same  as  two  copies 

of  the  roots  of  Eq  (90).  The  difference  is  that  the  A  coefficient  of  Eq  (90) 
is  equal  to  zero  because  the  eccentricity  of  satellite  2  is  equal  to  zero.  This 
reduces  Eq  (90)  to  a  1st  order  polynomial  (see  Eq  (94))  with  one  solution 
(see  Eq  (95)): 

B  cos(vJ  +  C  =  0  (94) 

cos(i/j)  =  -C/B  (95) 

However,  similar  to  above,  this  can  be  handled  by  checking  the  value  of  A 
before  solving  for  the  roots  of  Eq  (90).  If  A  equals  zero,  then  the  desired 
solutions  are  the  same  as  two  copies  of  the  single  root  of  Eq  (94)  (see  Eq 
(95)). 

The  procedures  above  yield  two  or  four  values  of  cos(i/})  that  are  the 
roots  of  Eq  (83),  or  Eq  (90)  when  the  absolute  value  of  yx  is  small  (on  the 

order  of  lO^km).  After  discarding  solutions  that  are  complex  or  have  an 
absolute  value  greater  than  1,  there  will  be  zero,  two,  or  four  valid  solutions 
left.  Ultimately,  when  there  are  two  valid  solutions,  the  close  approach 
boundaries  u2fn  and  must  be  found  that  meet  the  criteria  described  in 
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Eq  (18)  for  n  equal  to  one.  Sim&rily,  when  there  are  four  valid  solutions, 
the  close  approach  boundaries  i/J/n,  i/J/12,  v2f2V  and  v2f22  must  be  found  that 

meet  the  criteria  described  in  Eq  (18)  for  n  equal  to  two. 

Two  problems  remain.  The  first  problem  is  that  both  cos(i/2)  and 

cos(2rr-i/2)  are  equal  to  cos(i/2).  Given  that  0  =  cos“1[cos(i/J)|  ,  it  is  not 
known  whether  v2  =  0  or  i/2  =  27T-0.  Of  course,  v2  can  be  found 
through  the  following  procedure: 

if  d(i/j,0)  =  dra 

else 

V%  -  277-0 

Unfortunately,  this  method  requires  a  lot  of  CPU  time  to  implement,  so 
another  way  is  needed.  The  second  problem  is  that  once  the  values  of  v2 

are  found  where  d(uvv2)  =  dTO  ,  the  solutions  for  Vyiv  u2jiv  and 

i/2/22  are  still  not  known.  For  example,  with  two  v2  solutions,  there  is  no 
way  to  tell  which  of  the  two  solutions  is  i/J/n  or  v2fl2  without  some 
additional  work.  If  there  were  some  way  to  compute  i/2/n,  V2itv  an<* 

v2f22  directly  from  0,  then  both  problems  would  be  solved. 

Let  <t>x  through  equal  the  valid  solutions  for  cos(i/2),  such  that 

1  £  i  £  n 
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where 
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n  =  number  of  valid  solutions  for  cos(i/2) 

Let  0j  =  cos"1^)  where  1  ^  i  £  n  and  0  2*  0j  £  tt.  Note,  because  0j 
is  sorted  in  descending  order  (greater  i,  smaller  0j),  0j  will  be  sorted  in 
ascending  order  (greater  i,  larger  0j).  The  goal  is  now  is  to  find  some  way 
to  relate  0j  -  0B,  to  true  anomaly  of  the  close  approach  boundaries  of  each 
close  approach  region. 

This  process  is  simplified  considerably  by  redefining  close  approach  so 
that  a  close  approach  occurs  when  satellite  2  is  within  some  distance 
threshold  pdTH  of  the  projection  of  the  position  of  satellite  1  onto  the  orbital 

plane  of  satellite  2.  This  new  definition  of  close  approach  effectively  makes 
close  approach  a  two  dimensional  problem,  and  the  new  definition  of  close 
approach  is  completely  equivalent  to  the  old  definition,  so  long  as 

PdTH  =  (dTH  *  ZP'/5  (»«) 


where 

Zj  =  the  dbtance  from  satellite  1  to  the  orbital  plane  of  satellite  2 
(see  Eq  (69)) 

From  Eq  (69),  the  projection  of  the  position  vector  of  satellite  1  in  the 
perifocal  frame  of  satellite  2  on  to  the  orbital  plane  of  satellite  2  can  be 
represented  by 


ri/p?  =  < x.  yi  0 ) 


(97) 
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Two  final  terms  of  interest  are: 
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B,  =  {[^(l-ej-xj’+y*}*/*  (98) 

B2  =  {[■.(l+ej+xj’+y*}*  (99) 

where  Bt  is  the  distance  from  r  to  satellite  2’s  perigee,  and  B2  is  the 
distance  from  r2  to  satellite  2’s  apogee. 

Now  all  the  tools  are  in  place.  Given  the  results  of  Eqs  (96)  -  (99), 
there  are  three  basic  checks,  that  along  with  the  number  of  valid  solutions 
to  Eq  (83),  can  be  used  to  find  a  set  of  equations  relating  v^iv  v2f2V 

and  v2f22  to  0r  These  three  checks  are 

1.  Is  yj  >  0 

2.  Is  pdTH  >  B, 

3.  Is  pdTH  >  B2 

Since  the  results  of  each  check  is  either  true  or  false,  there  are  8  possible 
combinations  of  results.  Each  one  of  these  combinations  represents  a 
different  type  of  close  approach  which  requires  up  to  three  different  sets  of 
equations  to  represent  possible  cases  with  zero,  two,  and  four  valid  solutions. 
Table  1  lists  the  type  of  close  approach  that  corresponds  with  each  possible 
result  of  the  three  checks. 
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TABLE  1 
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Conditions  Required  for  Each  Type  of  Close  Approach 


Results  of  Close  Approach  Type  Checks 

Close  Approach 
Type 

y,  ^  o 

PdTH  *  B. 

PdTH  -  ^2 

1 

y,  *  o 

PdTH  *  B1 

PdTH  >  B» 

2 

y,  £  o 

PdTH  >  B1 

PdTH  *  B, 

3 

y,  ^  o 

PdTH  >  B1 

PdTH  >  B2 

4 

y,  >  o 

PdTH  ^  B1 

PdTH  *  Bj 

5 

y,  >  o 

PdTH  S  B1 

PdTH  >  BJ 

6 

y,  >  0 

PdTH  >  B1 

PdTH  *  B2 

7 

y  -  0 

PdTH  >  B1 

PdTH  >  B2 

8 

Appendix  A  contains  the  actual  equations  for  each  type  of  close 
approach.  As  a  general  convention  within  each  type  of  close  approach,  when 
there  are  two  valid  solutions,  both  v7f2l  and  are  set  to  zero.  When 

there  are  no  valid  solutions,  then  the  following  procedure  is  used  (also  see 
Eqs  (19)  and  (20)): 

(100) 
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if  pdfH  >  Bj  and  pd?H  >  B3 


V2fll  *“  ^ 


":/»  =  0 


l'vn  =  2r 


11P2  ~  ® 


and  pdTO  <  B, 


(101) 
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if  pdra  <  B, 
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Figure  2.  Example  of  a  Type  2  Close  Approach 

The  following  example  will  demonstrate  how  to  use  the  type  of  close 
approach  to  calculate  i/3/n,  i/?/ir  i/J/31,  and  from  04  through  0^ 

Inspection  of  Figure  2  reveals  that  yx  £  0,  pdTH  <  Br  and  pdTH  >  Br 

This  means  that  Figure  2  is  an  example  of  a  type  2  close  approach.  Since 
the  orbit  of  satellite  2  enters  and  exits  the  area  of  close  approach  once, 
there  are  two  valid  solutions.  Once  again  by  inspection,  0t  is  about  160°, 
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and  8j  is  about  170°.  For  two  valid  solutions,  Table  A-2  contains  the 
equations  to  calculate  vifu  and  from  0j  and  02: 


vVn  =  9s 


vm  =  2n  -  ei 


(102) 

(103) 


Applying  Eqs  (102)  and  (103),  i/J/u  equals  170°,  and  v2fli  equals  200°.  Since 
there  are  only  two  valid  solutions,  i/2/21  and  z/2/22  are  by  definition  equal  to 


Algorithm  Summary 

The  probability  of  close  approach  between  two  elliptical  orbits  can  be 
found  by  numerically  integrating  the  following  equation  (also  see  Eq  (22)): 


PCA  =  J  AM}(i/,)  (■ 


T1)  A 


(104) 


Where  for  each  numerical  integration  step,  dM^i/jJ/di/j  is  computed  using 
(also  see  Eqs  (25)  and  (42)) 


=  [(l-e^/d+e,)]1/4 


(105) 


dM(i/,)/di/l 


[cos1^,/2)  +  sin*(i/,/2)3* 


(106) 


and  AMjti/j)  is  computed  using  the  following  procedure: 


1.  Compute  the  positions  of  both  satellites  within  the  perifocal  frame 
of  satellite  2  using  Eqs  (43)  -  (62),  and  (68). 

2.  If  the  absolute  value  of  yx  is  small  (on  the  order  of  10“8  km)  and 
e}  is  non  zero,  then  find  the  roots  of  Eq  (90)  and  use  two  copies  of 
those  roots  to  obtain  four  solutions.  If  the  absolute  value  of  yt  is 
small  and  e2  is  zero,  then  find  the  root  of  Eq  (94)  and  use  two  copies 
of  that  root  to  obtain  two  solutions.  If  the  absolute  value  of  yx  is 
not  small,  find  the  roots  of  Eq  (83). 

3.  Discard  those  roots  that  are  complex,  or  those  with  absolute 
values  that  exceed  one.  The  roots  that  remain  are  valid  solutions. 

4.  Let  through  0^  equal  the  valid  solutions  for  cos(i/a)  such  that 

1  £  i  £  n 

where  n  is  the  number  of  valid  solutions.  Let  0j  =  cos“*(0j)  where 
1  £  i  £  n  and  0  &  Gj  £  7r. 

5.  Using  Eqs  (98)  and  (99),  compute  the  two  distance  bounds, 
and  Bj.  Use  Eq  (9G)  to  compute  the  close  approach  projected 
distance  threshold,  pdTH. 


6.  Perform  the  close  approach  checks,  and  locate  the  appropriate 

type  of  close  approach  in  Table  1.  Look  up  the  desired  type  of  close 
approach  in  Appendix  A,  then  using  the  number  of  valid  solutions, 
select  the  proper  equations  relating  i/J/n,  and  i/J/M  to  0r 

Compute  vvn,  and  i/„„. 

7.  Calculate  AM2(i/j)  with  the  following  equation: 

AM, (i/,)  =  M,(t/,/1})  -  M,(t/,/u)  +  M,(i/„„)  -  M,(i/}/„)  (107) 


where  (also  see  Eqs  (25)  and  (35)) 


0,  =  [(l-e,)/(l+e,)]‘/l 


(108) 


M,(y,)  =  2  tan(i/,/2)]  - 


e,  <?,  sin(t/,/2)  cos(t/,/2) 
(cos1(i/,/2)  -  fi\  sin‘(i/,/2)] 


(109) 
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IV.  Probability  of  Close  Approach  Between 
Satellites  in  Circular  Orbits 

Full  Circular  Orbit  Method 

The  purpose  of  this  section  is  to  develop  a  method  for  calculating  the 
probability  of  close  approach  between  two  satellites  in  circular  orbits. 

To  simplify  the  analysis,  the  close  approach  of  two  satellites  is 
redefined  to  be  whenever  the  angle  between  the  radius  vectors  of  satellite  1 
and  satellite  2  (angle  D)  is  less  than  or  equal  to  some  angle  threshold  DTH, 

where  Dra  is  equal  to  the  angle  between  the  radius  vector  of  satellite  1  and 

satellite  2  when  the  distance  between  the  two  satellites  is  equal  to  the 
distance  threshold  of  close  approach,  dra.  When  a  close  approach  is  possible 

(|r1— r2(  <  dTH),  the  plane  trigonometry  law  of  cosines  (7:196)  can  be  used  to 

solve  for  as  a  function  of  dTH: 

dTH  =  ri*  +  ri*  "  2r,r}cos(DTH) 
r  r*  +  r?  -  d*  ] 

Dth  —  cos  J  (110) 


where 


Tj  =  magnitude  of  the  radius  vector  of  satellite  1 
r}  =  magnitude  of  the  radius  vector  of  satellite  2 
dTH  =  distance  threshold  for  close  approach 
Dth  =  angular  distance  threshold  for  close  approach 

Note  that  using  an  angle  threshold  of  DTH  is  completely  equivalent  to  using 
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a  distance  threshold  of  dTO  so  long  as  both  satellites  are  in  circular  orbits 
and  Dth  is  computed  using  Eq  (110).  The  use  of  Dra  simplifies  the 

remaining  mathematics,  because  by  projecting  the  position  of  satellite  2  onto 
a  sphere  of  radius  rv  spherical  trigonometry  can  be  used  to  obtain  the  limits 

of  integration. 


Figure  3.  Spherical  Geometry  of  Circular  Orbit  Close  Approach 
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Figure  3  shows  the  orbital  path  of  satellite  1,  and  the  projection  of  the 
orbital  path  of  satellite  2.  For  circular  orbits,  the  position  of  zero  mean 
anomaly  is  arbitrary,  so  for  convenience,  the  mean  anomaly  of  both  satellites 
(Mj  and  M,)  are  assumed  to  be  zero  where  the  two  orbital  planes  cross  in 

the  northern  hemisphere.  When  the  mean  anomaly  of  satellite  1  is  known,  a 
close  approach  is  possible  whenever  some  portion  of  the  projection  of  satellite 
2’s  orbit  comes  within  DTH  (great  circle  arc)  of  satellite  1.  The  probability 

of  close  approach  (PCA)  can  be  determined  by  integrating  the  joint  density 
function  (see  Eq  (10))  over  the  region  of  Mj  and  M,,  where  D  is  less  than 
or  equal  to  DTH.  If  MJ2  and  Mn  are  the  unknown  limits  of  integration  over 
Mr  then  PCA,  expressed  in  terms  of  M12  and  Mu,  is  equal  to 
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AM,  =  M„  -  Mai 


(112) 


o  rM«  am. 

PCA  =  dM, 

i’  v/ 


(H3) 
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Two  more  things  are  needed  to  compute  PCA.  First,  an  equation  for 
AM,  as  a  function  of  Mj  is  needed.  Second,  the  integration  limits  over 
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must  be  found  for  which  a  close  approach  with  satellite  2  is  possible.  Using 
spherical  trigonometry  (7:198-200),  Eqs  (114)  and  (115)  can  be  derived  (see 
Figure  3): 


sinfMjJ/sinfrr^)  =  sin(x)/sin(8) 

sin(x)  =  sin(0)  sin(Mj)  (114) 

cos(Dra)  =  cos(x)cos(AM1/2) 

AM,  =  2  cos’^cosfD^J/costx)]  (115) 


where 

x  =  angle  between  the  radius  vector  of  satellite  1  and  the  orbital  plane 
of  satellite  2 

Given  the  trigonometric  relationship  (7:188) 

cos(x)  =  [l-sin^x)]1'1 

Eq  (114)  can  be  used  to  obtain  an  equation  for  cos(x): 

cos(x)  =  [1  -  sin1(e)sin,(Ml)]1^  (116) 

Substituting  Eq  (116)  into  Eq  (115),  yields  the  final  equation  for  AM,: 

AM,  =  2  cos•l{cos(DTH)/[l-sin,(e)*«l,(Ml))l',}  (117) 

The  integration  regions  can  be  found  through  Eq  (114),  by  replacing  x 
with  and  M(  with  M,  and  then  solving  for  M. 
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sin(DTH)  =  sin(O)  sin(M) 
sin(M)  =  sin(DTH)/sin(0) 

M  =  5in"l[sin(DTH)/sin(0)]  (118) 

If  8m(DTH)/sin(0)  is  less  than  1,  then  there  are  two  regions.  The  first 

region  is  from  -M  to  M,  and  the  second  is  from  7T-M  to  tt+M. 

If  sin(DTH)/sin(0)  is  greater  than  ly  then  the  integration  is  from  0  to 

2tr,  because  in  this  case,  for  0  <  M,  <  2tt  ,  there  is  always  some  chance  of 

a  close  approach  with  satellite  2. 

Assuming  that  there  are  some  places  in  satellite  l’s  orbit  where  there  is 
no  possibility  of  a  close  approach  by  satellite  2,  then  the  final  equation  for 

Pc*  “ 


There  is  no  closed  form  solution  for  the  equation  above,  so  numerical 
integration  must  be  used  to  obtain  the  final  solution  for  probability  of  close 
approach  between  two  satellites  in  circular  orbits. 
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Circular  Orbit  Approximation  Method 


In  the  previous  section  no  closed  form  solution  could  be  be  found  for 
PCA  between  two  satellites  in  circular  orbits.  However,  there  is  a  special 
case  which  does  have  a  closed  form  solution  for  PCA. 

Assume  that  satellite  1  and  satellite  2  are  in  circular  orbits  where  DTH 
is  small,  and  0  is  not  small.  This  also  implies  that  M2  and  Mj  (or  Mj-tt) 

are  also  small.  For  small  absolute  values  of  u,  the  following  approximations 
will  be  of  use  (with  3rd  order  effects  and  higher  discarded) (7:454-457): 

sin(u)  =  u 
cos(u)  =  1  -  u3/2 
(1  -  u3)1'2  =  1  -  u3/2 
1/(1  -  ur)  *  1  +  u3 

Eq  (117)  can  be  re-written  as 

cos(AM2/2)  =  cos(DTH)  /  [1  -  sin2(0)  sin3(M1)]1/3  (120) 

Substituting  the  small  value  approximations  into  Eq  (120)  and  then 
simplifying  yields 

(l-AMJ/8)  =  (l-D^/2)  [1  +  sin3(0)  Mj/2] 

1  -  AM3/8  -  1  4-  sin3(0)  Mj/2  -  D3H/2  -  D^inai{0^ 

^4th  Order 

AM{/8  =  D3h/2  -  sin3(0)  M*/2 
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AM,  =  2  [  D2ra  -  sin3(0)  U\  J1'3 
AM,  =  2  sin(O)  [  D3H/sin3(0)  -  MJ  ]1'3 


(121) 


When  Mj— 7T  is  very  small,  then  another  set  of  equations  must  be  used. 
If  Mj  =  7T-fu  ,  where  the  absolute  value  of  u  is  small,  then  the  following 
approximations  can  be  used: 

sin(Mj)  =  sin(7T-u) 

sin(Mj)  =  sin(7r)  cos(u)  -f  sin(u)  cos(7r) 
sin(Mj)  =  -sin(u) 
sin(Mj)  —  -u 

sin(Mj)  S  7T— Mj  (122) 

Substituting  the  small  value  approximations  (with  small  Mj-7t)  into  Eq 
(120),  results  in 

(l-AM3/8)  =  (l-DTH/2)  [1  4-  sin3(0)  (7T-Mj)3/2] 

1  -  AM»/8  =  1  +  sinJ(0)(Tr-Ml),/2  -  D^/2  - 

4th  Order 

AM*/8  =  Dsth/2  -  sin’(0)  (n-M,)72 
AM,  =  2  [  D’h  -  sins(0)  (n-M,)5  I1'5 

AM.  =  2  sin(0)  [  D*H/sins(0)  -  (rr-M,)1  l1*’  (123) 

By  substituting  the  small  value  approximations  into  Eq  (118),  the 

46 


limits  of  integration  can  be  found: 


M  =  DTH/sin(0) 


(124) 


Eqs  (121)  and  (122)  are  approximations  for  Eq  (117),  and  Eq  (124)  is  an 
approximation  for  Eq  (118).  Substituting  Eqs  (121),  (123),  and  (124)  into 
Eq  (119),  results  in: 


sin(0)  [  D|H/sinJ(0)  -  M]  ] W 

2P 


dMj 


p7T-fM 
*  rr-M 


sin(O)  [  DVsin*(0)  -  (tt-M,)’  ]*/* 
27? 


dMj 


(125) 


where  Eq  (124)  is  used  to  compute  M.  Given  that  (7:411) 

J  (a*-u*)I,a  du  =  [  u  (a*-u*)1^  4-  a*  sin^1(u/|a|)  ]/2  (126) 

Eq  (125)  can  be  directly  integrated  by  making  the  following  substitutions 
into  Eq  (126): 


For  small  Mji 
a  =  DTH/sin(0) 
u  =  M, 
du  =  dMj 


For  small  Mj-tr: 
a  =  Djg/suKO) 
u  =  tt-Mj 
du  =  -dMj 


Integrating  Eq  (125)  and  simplifying  results  in 
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47T*  L  1  \sm*(0) 


+ 


— k_  sin-i/ _ Mi _ \ 

sin’W  ViDn/BDOJi/J 


Dra/sin(e) 

-Dra/sin(e) 


+ 


/  7T— M,  V 

ain*(0)  Sm'\  |DTH/sin(e)|  j 


DL, 

UJ3L 


^+ptH/sin(9)] 

Tr-Pnj/sinCfl)] 


_  sin(6)  /  2tt  D’m  \ 
CA  47T1  \  sin,(9)  / 


D3 

^TH 

2 7T  sin(6) 


(127) 


Equ  (128)  is  the  closed  form  approximate  solution  for  P0A,  where  satellite  1 
and  satellite  2  are  in  circular  orbits,  DTH  is  small,  and  0  is  not  small. 
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To  verify  the  algorithms  from  the  previous  chapters,  three  computer 


programs  were  created. 

The  first  program  is  called  Statistical  Simulation  of  Probability  of  Close 
Approach  or  SSPCA.  SSPCA  queries  the  user  for  a  close  approach  distance 
threshold  and  for  a  set  of  orbital  parameters  for  satellite  1  and  satellite  2 
and  then  computes  the  probability  of  close  approach  between  both  satellites 
through  a  statistical  simulation.  The  first  step  in  this  process  is  to  select 
two  random  numbers  that  are  uniformly  distributed  between  0  and  2n  to 
represent  the  mean  anomalies  of  the  two  satellites.  For  the  selected  mean 
anomalies,  SSPCA  calculates  the  distance  between  the  two  satellites.  If  the 
computed  distance  between  the  two  satellites  is  less  than  or  equal  to  the 
input  distance  threshold,  then  a  close  approach  occurs.  This  process  is 
repeated  100,000  times,  and  a  count  is  kept  of  how  ma»;;  close  approaches 
occurred.  Eq  (128)  is  then  used  to  calculate  the  simulated  probability  of 
close  approach: 

Ps  =  nCA  /  100,000  (128) 

where  P$  is  the  simulated  probability  of  close  approach,  and  nCA  is  the 

number  of  close  approaches  that  occurred  in  the  simulation. 

The  second  program  is  called  Circular  Orbit  Probability  of  Close 
Approach  or  COPCA.  COPCA  first  queries  the  user  for  a  close  approach 


distance  threshold  and  for  a  set  of  orbital  parameters  for  two  circular  orbits. 
COPCA  then  computes  the  probability  of  close  approach  between  the  two 
satellites  using  both  the  full  and  the  approximate  circular  orbit  methods  of 
calculating  probability  of  close  approach  that  were  described  in  chapter  IV. 

The  third  program  is  called  Elliptical  Orbit  Probability  of  Close 
Approach  or  EOPCA.  EOPCA  queries  the  user  for  a  distance  threshold  and 
for  a  set  of  orbital  parameters  for  the  two  satellites  and  then  computes  the 
probability  of  close  approach  using  the  elliptical  orbit  method  of  calculating 
probability  of  close  approach  that  was  described  in  chapter  III. 

The  numerical  integration  in  both  the  COPCA  and  EOPCA  programs 
were  performed  using  Simpson’s  rule  with  an  step  size  of  approximately 
27t/10,000  radians. 

Statistical  Simulation  of  Probability  of  Close  Approach 

For  analysis  purposes,  assume  that  the  analytically  derived  probability 
of  close  approach  (PCA)  is  correct.  For  large  sample  sizes  and  values  of  both 

PCA  and  (1-PCA)  which  are  not  small,  the  number  of  close  approaches  that 

occur  in  the  statistical  simulation  can  be  approximated  by  a  normal 
distribution  with  a  mean  of  (5:225-226) 

nCA  “  ns  ^ca  (129) 


where 

nCA  -  mean  number  of  close  approaches  in  the  simulation. 
ns  =  number  of  samples  or  iteration?  in  the  simulation. 
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and  a  standard  deviation  of  (5:225-226) 

=  C  “s  PCA  0  -Pca)  ],/J  030) 

Likewise,  by  dividing  both  Eqs  (129)  and  (130)  by  ns,  Pg  can  be 
approximated  by  a  normal  distribution  with  a  mean  of  P0A,  and  a  standard 
deviation  of 


<7p  =  [  PCA  (1-P0A)  /  ns  ]V*  (131) 

Using  Eq  (131),  the  difference  between  Pg  (see  Eq  (128))  and  PCA  can 

now  be  found  in  terms  of  standard  deviations.  For  a  perfect  normal 
distribution,  the  absolute  value  of  the  difference  between  Pg  and  PCA  will  be 

less  than  .6745  a?  for  .5  of  the  simulation  runs,  and  less  than  1.96  ap  for 

.95  of  the  simulation  runs  (7:578).  These  two  thresholds  are  tests  that  can 
determine  how  well  the  simulated  solution  for  probability  of  close  approach 
matches  the  analytical  solutions  for  PCA  from  the  COPCA  and  EOPCA 

programs. 

Test  Cases 

Table  2  contains  a  list  of  the  orbital  parameters  of  satellite  1  and 
satellite  2  that  are  held  constant  through  all  test  cases.  Table  3  contains  a 
list  of  the  orbital  parameters  of  satellite  1  and  satellite  2  for  each  of  the  16 
test  cases  that  were  used  to  verify  the  probability  of  close  approach 
algorithms  described  in  chapters  III  and  IV.  Through  all  test  cases,  the 
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eccentricity  of  both  satellites  varied  from  0  to  .5,  with  the  eccentricity  of 
both  satellites  being  the  same  within  each  test  case.  For  convenience,  the 
longitude  of  the  ascending  node  of  both  satellites  and  the  inclination  of 
satellite  1  were  set  to  0°  so  that  the  inclination  of  satellite  2  would  equal  to 
the  angular  separation  of  the  two  orbital  planes. 


TABLE  2 

Orbital  Parameters  Held  Constant  Through  All  Test  Cases 


Orbital  Parameters 

Satellite  1 

Satellite  2 

Perigee  Radius 

7000  km 

7500  km 

Inclination 

0° 

varies 

Argument  of  Perigee 

0° 

90° 

Longitude  of  the  Ascending  Node 

0° 

0° 

Test  cases  1-8  (orbits  with  an  eccentricity  of  0  or  .1  )  were  run 
through  the  SSPCA  program  with  distance  thresholds  of  1000,  2000,  4000, 
8000,  12000,  and  20000  km.  Test  cases  9-16  (orbits  with  an  eccentricity  of 
.3  or  .5)  were  run  through  the  SSPCA  program  with  distance  thresholds  of 
4000,  8000,  12000,  and  2000U  km,  making  a  total  of  80  SSPCA  runs. 

All  circular  orbit  test  cases  (test  cases  1-4)  were  run  through  the 
COPCA  program  with  distance  thresholds  of  1000.  2000.  4000,  8000.  12000. 
20000  km.  making  a  total  of  24  COPCA  runs. 
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All  test  cases  were  run  through  the  EOPCA  program.  All  test  cases 
with  an  eccentricity  of  0  or  .1  (test  cases  1-8)  were  run  with  distance 
thresholds  of  1000,  2000,  4000,  8000,  12000,  and  20000  km.  All  test  cases 
with  an  eccentricity  of  .3  to  .5  (test  cases  9-16)  were  run  with  distance 
thresholds  of  4000,  8000,  12000,  and  20000  km,  making  a  total  of  80  EOPCA 
runs. 

Note  that  the  distance  thresholds  of  1000,  and  2000  km  were  only  used 
in  SSPCA,  COPCA,  and  EOPCA  runs  involving  test  cases  with  an 
eccentricity  of  .1  or  less.  This  is  because  when  the  distance  threshold  drops 
below  4000  km  and  the  test  case  eccentricity  is  .3  or  larger,  the  probability 
of  close  approach  is  generally  too  small  for  a  statistical  simulation  to  be  of 
much  value. 
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VI.  Results  and  Discussion 


Tables  B-l  and  B-2  in  Appendix  B  contain  the  test  results  of  the 
COPCA  program,  along  with  the  corresponding  simulation  results.  Table  4 
shows  how  well  the  simulated  PCA  matched  the  PCA  computed  in  COPCA 

for  all  runs  where  PCA  is  not  equal  to  1.  Runs  with  a  PCA  of  one  were 
excluded  from  Table  4  because  the  computed  PCA  always  equaled  the 
simulated  PCA  when  the  computed  PCA  was  equal  to  one,  and  because  <rp 
equals  zero  when  PCA  equals  one.  The  mean  error  listed  in  Table  4  is  the 
mean  of  the  error  (scaled  by  a  )  between  simulated  PCA  and  analytical  PCA 
for  all  COPCA  runs  with  a  PCA  less  than  one. 


TABLE  4 

Simulated  Probability  of  Close  Approach  Versus 
COPCA  Probability  of  Close  Approach 


Fraction  of 
COPCA  Runs 

With  Errors 

Normal 

Distribution 

Simulated 

^CA 

Less  Than  .6745  a 

.5000 

.5000 

Less  Than  1.96  a f 

.9500 

1.0000 

Mean  Error  (<7^) 

.0000 

.4445 

Tables  B  3  and  B  1  in  Appendix  B  contain  the  test  result  •  of  the 
' v'lJ'  '  -  “l“—  ■' —  - tach  is  comouted  usine 
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both  the  full  circular  orbit  method  and  the  circular  orbit  approximation 
method  described  in  chapter  IV.  For  all  runs  with  a  PCA  less  than  .01,  the 

circular  orbit  approximation  method  agreed  with  the  full  circular  orbit 
method  to  within  \% .  When  the  angle  between  the  orbital  planes  is  60°  or 
greater,  the  error  between  the  full  circular  orbit  method  and  the  circular 
orbit  approximation  method  was  generally  less  than  3%  when  the  computed 
PCA  was  less  than  5%.  As  expected  (see  Eq  (127)),  given  the  distance 

threshold,  the  error  in  the  approximation  method  is  inversely  proportional  to 
the  angle  between  the  two  orbital  planes.  Similarly,  given  the  angle  between 
the  orbital  planes,  the  error  in  the  approximation  method  is  directly 
proportional  to  the  distance  threshold. 


TABLE  5 

Simulated  Probability  of  Close  Approach  Versus 
EOPCA  Probability  of  Close  Approach 


Fraction  of 

EOPCA  Runs 

With  Errors 

Normal 

Distribution 

Simulated 

Pc* 

Less  Than  .6745  af 

.5000 

.5000 

Less  Than  1.96  a 

.9500 

.9657 

Mean  Error  (<7^) 

.0000 

.00"4 

Tables  B-5  through  B-10  in  Appendix  B  contain  the  test  results  of 
the  EOPCA  program,  along  with  the  corresponding  simulation  results.  Table 
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5  shows  how  well  the  simulated  PCA  matched  the  PCA  computed  in  EOPCA 
for  all  runs  where  PCA  is  not  equal  to  1,  As  with  the  COPCA  test  results, 
runs  with  a  PCA  of  one  were  excluded  from  Table  5,  because  the  computed 
PCA  always  equaled  the  simulated  PCA  when  the  computed  PCA  was  equal  to 
one,  and  because  a*  =  0  when  PCA  =  1  .  The  mean  error  listed  in 
Table  5  is  the  mean  of  the  error  (scaled  by  aj  between  simulated  PCA  and 
analytical  PCA  for  all  EOPCA  runs  with  a  PCA  less  than  one. 

The  probability  of  close  approach  computed  by  both  the  COPCA  and 
EOPCA  programs  favorably  matches  the  simulated  probability  of  close 
approach  computed  by  the  SSPCA  program.  However,  the  mean  error  for 
the  20  COPCA  runs  with  a  computed  PCA  less  than  one  (see  Table  4) 

indicates  that  a  possible  bias  exists  between  simulated  PCA  and  the  analytical 
PCA  computed  by  COPCA.  Similar  biases  exist  within  the  EOPCA  test 

results  when  EOPCA  runs  with  only  the  same  eccentricity  are  examined. 
When  all  72  EOPCA  runs  with  a  computed  PCA  less  than  one  (see  Table  5) 

are  considered,  there  does  not  appear  a  bias.  Given  the  limited  number  of 
circular  orbit  test  cases,  the  small  bias  in  the  COPCA  test  results  is  not 
considered  significant. 
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VII.  Suggestions  and  Recommendations 


To  develop  a  method  to  calculate  the  probability  of  close  approach 
between  two  satellites,  two  major  assumptions  were  made.  First,  all  orbital 
elements,  except  for  true  (or  mean)  anomaly,  were  assumed  to  be  constant 
over  time.  Second,  the  mean  anomalies  of  both  satellites  were  assumed  to 
be  independent  random  variables  that  are  uniformly  distributed  between  0 
and  2tt.  While  these  assumptions  greatly  simplified  the  derivation  of 
probability  of  close  approach,  they  also  limited  its  usefulness. 

Three  follow'- up  studies  are  recommended.  The  goal  of  the  first 
study  would  be  to  develop  methods  to  calculate  probability  of  close  approach 
between  two  satellites  in  elliptical  orbits,  where  there  are  linear  perturbations 
to  the  argument  of  perigee  and  the  longitude  of  the  ascending  node  of  both 
satellites.  In  this  case,  both  argument  of  perigee  and  longitude  of  the 
ascending  node  would  be  treated  as  linear  functions  of  time.  The  goal  of 
the  second  study  would  be  to  develop  methods  to  calculate  probability  of 
close  approach  between  two  satellites  in  elliptical  orbits  where  the  duration 
of  the  specified  time  interval  is  too  short,  or  the  ratio  of  the  orbital  periods 
of  the  two  satellites  are  such  that  the  mean  anomalies  of  both  satellites  are 
not  independent.  The  goal  of  the  final  study  would  be  to  find  ways  to 
reduce  the  computational  expense  involved  in  calculating  probability  of  close 
approach.  In  this  paper,  numerical  integration  was  used  to  directly  calculate 
the  probability  of  close  approach  between  two  satellites.  This  approach  can 
be  used  to  accurately  calculate  the  probability  of  close  approach  between  any 
two  satellite:  with  eccentricities  leos  than  1.0.  but  it  can  be  computational!) 
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expensive.  For  two  satellites  with  an  eccentricity  less  than  .3,  it  is  possible 
that  a  series  approximation  for  probability  of  close  approach  could  be  found 
that  would  offer  acceptable  precision,  and  at  much  less  computational 
expense  than  methods  that  use  numerical  integration.  For  two  satellites  in 
circular  orbits,  even  simpler  series  approximations  for  probability  of  close 
approach  could  ue  possible.  Both  forms  of  series  approximations  for 
probability  of  close  approach  merit  further  investigation. 
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Type  1  Close  Approach 
A  type  1  close  approach  occurs  when 


1.  y,  £  0 
2-  pdTH  £  Bj 

3.  pdra  £  B, 


When  there  are  0  valid  solutions,  then 


o 

II 

V* 

^2/21  ~ 

o 

II 

« 

M 

vyn  = 

When  there 

are  2  valid  solutions,  then 

Vvu  =  2n  ~  es 

^2/21  = 

vun  =  2^-0, 

^2/22  = 

When  there 

are  4  valid  solutions,  then 

"*/..  =  0S 

VH%1  = 

vU\i  =  ®i 

^2/22  = 
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Type  2  Close  Approach 
A  type  2  dose  approach  occurs  when 


1.  y,  S  0 

2.  pdTH  £  B, 

3-  PdTH  >  BJ 

It  is  not  possible  to  have  0  valid  solutions  in  a  type  2  close  approach. 
When  there  are  2  valid  solutions,  then 


2/11 

^2/21 

3/12 

=  2tt  -  e, 

V2/22 

When  there  are  4  valid  solutions,  then 


'■in 

=  », 

‘'i/n  =  2lt  - 

-0. 
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Type  3  Close  Approach 
A  type  3  close  approach  occurs  when 


1.  y,  ^  0 

2.  pdTH  >  B, 
3-  pdfj,  £  B, 


It  is  not  possible  to  have  0  valid  solutions  in  a  type  3  close  approach. 
When  there  are  2  valid  solutions,  then 


Vl,U 

=  -0, 

vmx 

=  0 

=  0, 

vxin 

=  0 

valid  solutions,  then 

vVu 

=  -0, 

Vvn 

II 

CP 

vV\t 

=  e, 

viir> 

=  0, 
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Type  4  Close  Approach 
A  type  4  close  approach  occurs  when 


i 

i 


2 

V 

v 

H 

I 

! 


V 


V 


v 


1.  y,  £  0 

2.  pdTH  >  B, 

3.  pdTH  >  Bs 


When  there  are  0  valid  solutions,  then 


vV\\  ~  ® 

^2/21 

II 

^2/22 

When  there 

are  2  valid  solutions,  then 

vtiu  =  e»  -  2n 

^2/21 

~  ®i 

V2f22 

When  there 

are  4  valid  solutions,  then 

"vij  =  0,  -  2rr 

vw 

vt/u  =  ®i 

V2f22 
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Type  5  Close  Approach 
A  type  5  close  approach  occurs  when 


1.  y,  >  o 

2.  pdTH  S  Bj 

3.  pdra  *  Bs 


When  there  are  0  valid  solutions,  then 


"yn 

=  0 

V2/2l 

=  0 

vt/n 

=  0 

V2f22 

=  0 

valid  solutions,  then 

vtin 

=  ®, 

V2!2i 

=  0 

vVn 

c* 

O 

ii 

V2!22 

=  0 

valid  solutions,  then 

vWl 

=  e, 

V2f2\ 

=  2rr  ~  0f 

vtin 

=  ®« 

V2t22 

=  27T  -  03 
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Type  6  Close  Approach 
A  type  6  close  approach  occurs  when 


1.  y,  >  0 

2-  pdTH  ^  B, 
3.  pdra  >  B, 


It  is  not  possible  to  have  0  valid  solutions  in  a  type  6  close  approach. 


When  there  are  2  valid  solutions,  then 


vwi  =  e. 


Vvn  =  2?r  -  0» 


When  there  are  4  valid  solutions,  then 


^2/n 


Vva  =  *  -  °4 
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vVn  =  0 


"t/M  =  0 


"a/a.  =  27r  *  8. 


^aa  =  2n  -  ®a 
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Type  7  Close  Approach 
A  type  7  close  approach  occurs  when 


1.  yx  >  0 

2.  pdra  >  B, 

3.  pdTH  £  B, 


It  is  not  possible  to  have  0  valid  solutions  in  a  type  7  close  approach. 
When  there  are  2  valid  solutions,  then 


^2/11 

=  -  », 

Vt/U 

=  0 

Vll\2 

=  6: 

^2/22 

=  0 

valid  solutions,  then 

vy\\ 

= 

^2/21 

=  2tt  - 

-  », 

^2/12 

^2/22 

=  2tt  - 

-  o, 

k  i .  29  ft 


I 


A  type  8  close  approach  occurs  when 


1.  y,  >  0 

2.  pdra  >  5, 

3.  pdTH  >  B, 


When  there  are  0  valid  solutions,  then 


o 

II 

*'*/*! 

=  0 

vm..m  =  2rr 

vm  M 

=  0 

2/12 

1/22 

When 

there 

are  2 

valid  solutions,  then 

M 

II 

1 

JO 

^2/2 1 

=  0 

=  27r  "  e* 

^2/21 

0 

When 

there 

are  4 

valid  solutions,  then 

II 

t 

CD 

=  2rr 

l/ll  1 

2/21 

*'i/n  =  2r:  -  e4 

"l/li 

=  2n 
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TABLE  B-l 


COPCA  Test  Results  For  Distance  Thresholds  Less  Than  4000  km 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

Pea 

Simulation 

Simulated 

Pc* 

Error  (af) 

Error  (%) 

1000 

0° 

.038068 

.000605 

.038510 

.7306 

1.16 

2000 

0° 

.085327 

.000883 

.086210 

1.0000 

1.03 

1000 

30° 

.004580 

.000214 

.004670 

.4206 

1.97 

2000 

30° 

.023638 

.000480 

.024010 

.7750 

1,57 

1000 

60° 

.002632 

.000162 

.002790 

.9753 

6.00 

2000 

60° 

.013287 

.000362 

.013250 

-.1022 

-.28 

1000 

90° 

.002278 

.000151 

.002370 

.6093 

4.04 

2000 

90° 

.011471 

.000337 

.011100 

-1.1009  i 

-3.23 

TABLE  B-2 


COPCA  Test  Results  For  Distance  Thresholds  Not  Less  Than  4000  km 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

Pca 

Simulation 

Simulated 

PCA 

Error  (o^) 

Error  (%) 

4000 

0° 

.176602 

.001206 

.178280 

1.3914 

.95 

8000 

0° 

.371487 

.001528 

.373250 

1.1533 

.47 

12000 

0° 

.620316 

.001535 

.620220 

-.0625 

-  02 

20000 

0° 

1.000000 

.000000 

1.000000 

— 

.00 

4000 

30° 

.128739 

.001059 

.129460 

.6808 

.56 

8000 

30° 

.361425 

.001519 

.363150 

1.1356 

.48 

12000 

30° 

.629654 

.001527 

.629590 

-.0419 

-.01 

20000 

30° 

1.000000 

.000000 

1.000000 

- 

.00 

4000 

60° 

.058214 

.000740 

.058130 

-.1135 

-.14 

8000 

60° 

.315640 

.001470 

.316940 

.8840 

.41 

12000 

60° 

.671401 

.001485 

.671640 

.1609 

.04 

20000 

60° 

1.000000 

.000000 

1.000000 

- 

.00 

4000 

90° 

.049658 

.000687 

.050050 

.5691 

.79 

8000 

90° 

.233877 

.001339 

.233390 

-.3637 

-.21 

12000 

90° 

.754522 

.001361 

.755730 

.8876 

.16 

20000 

90° 

1.000000 

.000000 

1.000000 

- 

.00 

TABLE  B-3 


COPCA  Test  Results  Using  Circular  Orbit  Approximation 
For  Distance  Thresholds  Less  Than  4000  km 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

Pci 

COPCA 

Approximate 

Pci 

Error  (%) 

1000 

30° 

.004580 

.004553 

-.59 

2000 

30° 

.023638 

.022873 

-3.24 

1000 

60° 

.002632 

.002629 

-.11 

2000 

60° 

.013287 

.013206 

-.61 

1000 

90° 

.002278 

.002276 

-.09 

2000 

90° 

.011471 

.011437 

-.30 
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TABLE  B-4 


COPCA  Test  Results  Using  Circular  Orbit  Approximation 
For  Distance  Thresholds  Not  Less  Than  4000  km 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

Pca 

COPCA 

Approximate 

Pca 

Error  (%) 

4000 

30° 

.128739 

.097981 

-23.89 

8000 

30° 

.361425 

.433547 

19.95 

12000 

30° 

.629654 

- 

— 

20000 

30° 

1.000000 

- 

- 

4000 

60° 

.058214 

.056569 

-2.83 

8000 

60° 

.315640 

.250309 

-20.70 

12000 

60° 

.671401 

.697936 

3.95 

20000 

60° 

1.000000 

— 

— 

4000 

90° 

.049658 

.048990 

-1.35 

8000 

90° 

.233877 

.216774 

-7.31 

12000 

90° 

.754522 

.604430 

-19.89 

20000 

90° 

1.000000 

- 

- 

TABLE  B-5 


EOPCA  Test  Results  For  Circular  Orbits  With 
Distance  Thresholds  Less  Than  4000  km 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

^CA 

Simulation 

Sim  u  nited 

^CA 

Error  (crp) 

Error  {%) 

1000 

0° 

.038i  J8 

.000605 

.038510 

.7306 

1.16 

2000 

0° 

.085327 

.000883 

.086210 

1.0000 

1.03 

1000 

30° 

.004580 

.000214 

.004670 

.4206 

1.97 

2000 

30° 

.023638 

.000480 

.024010 

.7750 

1.57 

1000 

60° 

.002632 

.000162 

.002790 

.9753 

6.00 

2000 

60° 

.013287 

.000362 

.013250 

-.1022 

-.28 

1000 

90° 

.002278 

.000151 

.002370 

.6093 

4.04 

2000 

90° 

.011471 

.000337 

.011100 

-1.1009 

-3.23 
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TABLE  B-6 


EOPCA  Test  Results  For  Circular  Orbits  With 
Distance  Thresholds  Not  Less  Than  4000  km 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

PCA 

Simulation 

Simulated 

PCA 

Error  (crj 

Error  (%) 

4000 

0° 

.176602 

.001206 

.178280 

1.3914 

.95 

8000 

0° 

.371487 

.001528 

.373250 

1.1538 

.47 

12000 

0° 

.620316 

.001535 

.620220 

-.0625 

-.02 

20000 

0° 

1.000000 

.000000 

1.000000 

- 

.00 

4000 

30° 

.128739 

.001059 

.129460 

.6808 

.56 

8000 

30° 

.361425 

.001519 

.363150 

1.1356 

.48 

12000 

30° 

.629654 

.001527 

.629590 

-.0419 

-.01 

20000 

30° 

1.000000 

.000000 

1.000000 

— 

.00 

4000 

60° 

.058214 

.000740 

.058130 

-.1135 

-.14 

8000 

60° 

.315640 

.001470 

.316940 

.8840 

.41 

12000 

60° 

.671401 

.001485 

.671640 

.1609 

.04 

20000 

60° 

1.000000 

.000000 

1.000000 

- 

.00 

4000 

90° 

.049659 

.000687 

.050050 

.5691 

.79 

8000 

90° 

.233877 

.001339 

.233390 

-.3637 

“.21 

12000 

90° 

.754522 

.001361 

.755730 

.8876 

.16 

20000 

90° 

1.000000 

.000000 

1.000000 

- 

.00 

TC 

I  o 


TABLE  B-7 


EOPCA  Test  Results  For  Eccentricity  of  .1  With 
Distance  Thresholds  Less  Than  4000  km 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

^CA 

Simulation 

Simulated 

PCA 

Error  (<Tp) 

Error  (%) 

1000 

0° 

.021034 

.000454 

.021520 

1.0705 

2.31 

2000 

0° 

.067767 

.000795 

.068180 

.5195 

.61 

1000 

30° 

.002449 

.000156 

.002030 

-2.6859 

-17.11 

2000 

30° 

.016417 

.000402 

.016320 

-.2413 

-.59 

1000 

60° 

.001414 

.000119 

.001170 

-2.0504 

-17.26 

2000 

60° 

.009292 

.000303 

.008960 

-1.0957 

-3.57 

1000 

90° 

.001225 

.000111 

.000950 

-2.4775 

-22.45 

2000 

90° 

.008027 

.000282 

.007530 

-1.7624 

-6.19 
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TABLE  B-8 


EOPCA  Test  Results  For  Eccentricity  of  .1  With 
Distance  Thresholds  Not  Less  Than  4000  km 


Distance 

Threshold 

(km) 

Angie 

Between 

Orbital 

Planes 

COPCA 

^CA 

Simulation 

Simulated 

^CA 

Error  (crf) 

Error  (%) 

4000 

0° 

.154566 

.001143 

.155130 

.4934 

.36 

8000 

0° 

.327955 

.001485 

.329990 

1.3704 

.62 

12000 

0° 

.533222 

.001578 

.533210 

-.0076 

-.00 

20000 

0° 

1.000000 

.000000 

1.000000 

- 

.00 

4000 

30° 

.091240 

.00091i 

.091410 

.1866 

.19 

8000 

30° 

.313524 

.001467 

.315190 

1.1357 

.53 

12000 

30° 

.535690 

.001577 

.535530 

-.1015 

-.03 

20000 

30° 

1.000000 

.000000 

1.000000 

- 

.00 

4000 

60° 

.045010 

.000656 

.045230 

.3354 

.49 

8000 

60° 

.234065 

.001339 

.233720 

-.2577 

-.15 

12000 

60° 

.546098 

.001574 

.547590 

.9479 

.27 

20000 

60° 

1.000000 

.000000 

1.000000 

i 

- 

.00 

4000 

90° 

.038555 

.000609 

.038810 

.4187 

.66 

8000 

90° 

.179613 

.001214 

.178870 

-.6120 

-.41 

12000 

90° 

.581394 

.001560 

.582690 

.8308 

no 

20000 

90° 

1.000000 

.000000 

1.000000 

- 

.ou 
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TABLE  B-9 


EOPCA  Test  Results  For  Eccentricity  of  .3 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

PCA 

Simulation 

Simulated 

P0A 

Error  (<rp) 

Error  {%) 

4000 

0° 

,067156 

.000791 

.066410 

-.9431 

-1.11 

8000 

0° 

.226341 

.001323 

.226230 

-.0839 

-.05 

12000 

0° 

.374049 

.001530 

.375610 

1.0203 

.42 

20000 

0° 

.762220 

.001346 

.762770 

.4086 

.07 

4000 

30° 

.027733 

.000519 

.027410 

-.6224 

-1.16 

8000 

30° 

.192271 

.001246 

.191640 

-.5064 

-.33 

12000 

30° 

.363538 

.001521 

.364400 

.5667 

.24 

20000 

30° 

.773937 

.001323 

.774410 

.3575 

.06 

4000 

60° 

.013120 

.000360 

.013530 

1.1389 

3.13 

8000 

60° 

.109096 

.000986 

.108750 

-.3509 

-.32 

12000 

60° 

.302533 

.001453 

.301440 

-.7522 

-.36 

20000 

60° 

.819091 

.001217 

.820710 

1.3303 

1.98 

4000 

90° 

.010933 

.000329 

.011360 

1.2979 

3.91 

8000 

90° 

.090037 

.000905 

.090000 

-.0409 

-.04 

12000 

90° 

.269083 

.001402 

.267690 

-.9936 

-.52 

20000 

90° 

.849413 

.001131 

.851050 

1.4474 

1.93 
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TABLE  B-10 


EOPCA  Test  Results  For  Eccentricity  of  .5 


Distance 

Threshold 

(km) 

Angle 

Between 

Orbital 

Planes 

COPCA 

^CA 

Simulation 

Simulated 

^CA 

Error  (<Tp) 

Error  (%) 

4000 

0° 

.015125 

.000386 

.014750 

-.9715 

-2.48 

8000 

0° 

.066378 

.000787 

.065430 

-1.2046 

-1.43 

12000 

0° 

.185335 

.001229 

.185330 

-.0041 

-.00 

20000 

0° 

.434822 

.001568 

.435010 

.1199 

.04 

4000 

30° 

.002293 

.000151 

.002270 

-.1523 

-1.00 

8000 

30° 

.042547 

.000638 

.041490 

-1.6567 

-2.48 

12000 

30° 

.150654 

.001131 

.150520 

-.1185 

-.09 

20000 

30° 

.427634 

.001564 

.427660 

.0166 

.01 

4000 

60° 

.000682 

.000083 

.000610 

-.8675 

-10.56 

8000 

60° 

.010386 

.000321 

.010230 

-.4860 

-1.50 

12000 

60° 

.108902 

.000985 

.108390 

-.5198 

-.47 

20000 

60° 

.390383 

.001543 

.389900 

-.3130 

-.12 

4000 

90° 

.000434 

.000066 

.000430 

-.0606 

-.92 

8000 

90° 

.008820 

.000296 

.008590 

-.7770 

-2.61 

12000 

90° 

.085002 

.000882 

.085920 

1.0408 

1.08 

20000 

90° 

.341475 

.001500 

.338930 

-1.6967 

-.75 
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Program  Listings 
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oooooo  o  oooooo 


. c 

c 

SSPCA  -  Statistical  Simulation  of  the  Probability  C 

of  Close  Approach  C 

C 

. - . C 

INTERFACE  TO  SUBROUTINE  RNDINICC3 
END 

INTERFACE  TO  REAL*8  FUNCTION  RANDOM [Cl 
END 


PROGRAM  SSPCA 

CHARACTER  YESNO*4 

LOGICAL*4  DEBUG.  OBSCUR 

INTEGERS  I.  J.  K.  ONE.  T NO.  CA.  OCA 

REAL*8  DX.  DY.  DZ.  DR 

REAL*8  Ml.  Al.  El.  INC1 .  NU1.  ARGPA1 .  LONAN1 .  XI (4) 

REAL*8  M2.  A2.  E2.  INC2.  NU2.  ARGPA2 .  L0NAN2.  X2(4) 

REAL*8  DTH,  RANDOM,  RANGE.  DBLE 

REAt*8  HALFPI . PI . TNOPI . DEGRAD . RADDEG 

REAL*8  ER.  DU.  TU 

COMMON  /ADMIN/  DEBUG 

COMMON  /ASTRO/  ER.  DU.  TU 

COMMON  /COMPAR/  DX.  DY.  DZ.  DR 

COMMON  /CONST/  HALFPI.  PI.  TWOPI .  DEGRAD,  RADDEG 

COMMON  /SAT17  Ml.  Al .  El.  INC1 .  NU1 .  ARGPA1 .  LONAN1 .  XI 

COMMON  /SAT2/  M2.  A2.  E2.  INC2.  NU2.  ARGPA2.  LONAN2 .  X2 

PARAMETER  (ONE-1 .  TNO-2) 

ER  -  6378 . 145D0 
DU  -  6378.145D0 
TU  •  806.81 18744D0 
P I -3 . 1 4 1 5926535897 93D0 
HALFPI-PI/2.0D0 
TNOPI-2.0D0*PI 
DEGRAD-PI/180. 0D0 
RADDEG- 1 80 . ODO/P I 


* — — - — - — - * 

•  Initialize  Random  Number  Generator  * 

*  - - - — - - - * 

CALL  RNDINI 

* - * 

*  Input  Orbital  Elements  of  Sat  1  * 

*  - - * 

«RITE(*.1000)  ONE 

READ  (•.1010)  Al 
NRITEC* . 1020)  ONE 
READ  '♦.1010)  El 


8! 


OOOO  OOOO  O  O  O  o  o  o 


WRITE (♦.1030)  ONE 
READ  (♦.1010)  INC1 
WRITE (♦. 1040)  ONE 
READ  (♦.1010)  ARGPA1 
WRITE (♦.1050)  ONE 
READ  (♦,1010)  LONAN1 


* - * 

♦  Input  Orbital  Elements  of  Sat  2  ♦ 

* — — - — — — * 


WRITER.  1000)  TWO 
READ  (♦.1010)  A2 
WRITE (♦.1020)  TWO 
READ  (♦.lOlO)  E2 
WRITE(* . 1030)  TWO 
READ  (♦.1010)  INC2 
WRITER  .1040)  TWO 
READ  (♦.1010)  ARGPA2 
WRITE (♦. 1050)  TWO 
READ  (♦.1010)  L0NAN2 


* — - - — — * 

♦  Input  Simulation  Limits  ♦ 

*  - - — - * 


WRITER.  1060) 

READ  (♦ . 1010)  DTH 
WRITE (♦.1070) 

READ  (♦ . 1080)  ITER 
WRITER.  1090) 

READ  (♦.1100)  YESNO 

IF  ( YESNO (1 : 1) .EG.  *Y * .OR. YESNO (1 : 1) .EQ. *y ')  THEN 
DEBUG  •  .TRUE. 

ELSE 

DEBUG  -  .FALSE. 

ENDIF 

* — - - - - * 

♦  Convert  the  input  angles  from  deg  ♦ 

•  to  rad  ♦ 

* — — — — - - — * 

I NCI  •  INC1  ♦  DEGRAD 
ARGPA1  •  ARGPAl  •  DEGRAD 
LONAN1  -  LON AN 1  ♦  DEGRAD 
INC2  «  INC2  ♦  DEGRAD 
ARGPA2  «  ARGPA2  ♦  DEGRAD 
L0NAN2  -  LONAN2  *  DEGRAD 


• - - ♦ 

♦  Clear  the  two  event  counters:  ♦ 

♦  OCA  «  Obscured  Close  Approach  • 

♦  CA  *  Close  Approach  * 


82 


* 

* 


c 


c 

c 

c 


100 


c 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 

1080 

1090 

1100 

1110 

1120 


♦  — 

OCA 

CA 


0 

0 


* - 

*  Start  Simulation 


♦ 


* — - — - — - * 

DO  100  1*1. ITER 
Ml  *  TWOPI  *  RANDOM () 

M2  -  TWOPI  ♦  RANDOM () 

DR  -  RANGE () 

IF  (DEBUG)  WRITE  (*.1110)  Ml,  M2.  DR 
IF  (DR  .LE.  DTH)  THEN 
IF  (QBSCURO)  THEN 
OCA  *  OCA  ♦  1 
ELSE 

CA  -  CA  *  1 
ENDIF 
ENDIF 
CONTINUE 

WRITE(* , 1120)  Al.  A2,  El.  E2.  INC1*RADDEG.  INC2*RADDEG , 


1  ARGPA1*RADDEG,  ARGPA2*RADDEG .  L0NAN1*RADDEG . 

2  L0NAN2*RADDEG .  DTH.  ITER. 

3  OCA.  DBLE(OCA)/DBLE(ITER) . 

4  CA.  DBLE (CA) /DBLE (ITER) 

* - - — - - - — + 


FORMAT  (/*  Incut  the  semi-major  axis  of  Sat  Ml. 

1  ’  in*  km:  ’\) 

FORMAT  (FI 6. 12) 

FORMAT  (’  Input  the  eccentricity  of  Sat  Ml.*:  *\) 

FORMAT  (*  Input  the  incllnition  of  Sat  Ml. 

1  ’  in  degrees:  *\) 

FORMAT  (*  Input  the  argument  of  perapsis  of  Sat  Ml. 

1  *  in  degrees:  *\) 

FORMAT  ('  Input  the  longitude  of  the  ascending  node  *. 

1  *of  Sat  Ml.1  in  degrees:  *\) 

FORMAT  (/*  Input  distance  threshold  in  km:  *\) 

FORMAT  (*  Input  the  desired  number  of  iterations  \ 

1  * (7  digits  max) :  *\) 

FORMAT  (17) 

FORMAT  (*  Run  SSPCA  in  DEBUC  mode  (Y/N)?  *\) 

FORMAT  (A4) 

FORMAT  (/’  Ml*’ ,F12. 10.4X.  *M2*’ . F12. 10.4X. * RANGE*' .F12.6) 
FORMAT  (/T21 .  ‘Sat  l*.T4l.*Sat  2*/Tl6. * . *. 

1  T36 .  * . *//’  a  (km)  *  .TIC .F16 , 10 ,4X . 

2  F16.10//*  e' .T16.F16. 10.4X.F16. 10// 

3  *  inc  (deg) ’.T16.F16.10.4X.F16. 10// 

4  *  argument  of ' .T16.F16. 10.4X.F16. 10/'  perigee  (deg)*// 
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5  ’  long  of  asc*,T16,F16.10.4X.F16.10/’  node  (deg)’/// 

6  ’  Distance  Threshold  for  Close  Approach:  * . 

7  F16 .  10 ,  *  km’/// 

6  T48, ’Fraction’/’  Iterations  ’.I7.T46,’ - *// 

7  ’  Close  ’.T14. 17,  ’—Obscured  by  the  Earth— ’  .F12. 10/ 


8  ’  Approaches  */T14 . 17,  * - Unobscured 

9  F12.10//) 

END 

* - * 

♦  * 

*  Function  OBSCURO  * 

*  * 

* - * 

LOG  I  CAL*  4  FUNCTION  OBSCURO 

C  * . . . . — . * 

LOGICAL*4  DEBUG 
INTEGER*4  I.  J.  K 


REAL*8  NDRANG.  DX.  DY.  DZ.  DR,  ODTH,  RADICAL 
REAL*8  Ml .  A1 .  El ,  INC1 .  NU1 .  ARGPA1 .  LONAN1 .  XI (4) 

REAL*8  M2.  A2.  E2.  INC2,  NU2.  ARGPA2.  L0NAN2,  X2C4) 

REAL*8  HALFPI . PI .TWOPI . DEGRAD . RADDEG 

REAL*9  ER.  DU.  TU 

COMMON  /ADMIN/  DEBUG 

COMMON  /ASTRO/  ER.  DU.  TU 

COMMON  /COMPAR/  DX.  DY.  DZ.  DR 

COMMON  /CONST/  HALFPI.  PI.  TWOPI.  DEGRAD.  RADDEG 

COMMON  /SAT  1 /  Ml.  A1 .  El.  INCl ,  NU1 .  ARGPA1 .  LONAN1 .  XI 

COMMON  /SAT2/  M2.  A2.  E2.  INC2.  NU2.  ARGPA2.  L0NAN2.  X2 

C  * - - - - — * 

NDRANG  -  -IDO  *  (  X1(1)*DX  ♦  IM'Z'ui  ♦  X1(3)*DZ  ) 

NDRANG  -  NDRANG  /  XI  (4)  '  DR 

IF  (DABSCNDRANG)  .GT.  IDO)  NDRANG  »  DSIGNC  IDO.  NDRANG) 
NDRANG  -  DACGS(NDkANC) 

IF  (DEBUG)  NR! TEC*. 1000)  NDRANG*RADDEG 
IF  (NDRANG. GT. HALFPI  OR.  XI (4) *DSIN (NDRANG) .GE.ER)  THEN 
OBSCUR  •  .FALSE. 

ELSE 

ODTH  "  XI (4)  *  DCOS(NDRANG) 

RADI  L  «  ODTH  •  ODTH  -  X1(4)*X1(4)  *  ER*ER 
IF  (RADICAL  .GE.  0D0)  THEN 
ODTH  -  ODTH  -  DSQRT (RADICAL) 

ELSE 

IF  (  DSQRT (  DABS (RADICAL)  )  .GT.  5D-11  )  THEN 
WRITE  (• . 10!0>  ODTH. 

1  -IDO  *  DSQRT (-RADICAL) 

EKDIF 

ENDIF 
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IF  (DEBUG)  WRITE  (*,1020)  DR,  ODTH. 

1  DSIGN (DSQRT (DABS (RADICAL) ) , RADICAL) 

IF  (  DR  .GT.  ODTH)  THEN 
OBSCUR  -  .TRUE. 

ELSE 

OBSCUR  *  .FALSE. 

ENDIF 

ENDIF 

RETURN 

1000  FORMAT  (//’  The  Nadir  Angle  from  Sat  1  to  Sat  2  \ 

1  ’is  (deg) : * ,F16. 10) 

1010  FORMAT  (//’  Error!  Negative  Radical.  ODTH*’ ,F16. 10, 

2  4X,  ’RADICAL*’,  1P.D12.5.0P) 

1020  FORMAT  (//’  DR-’ ,F16. 10. 4X. ’ODTH*’ ,F16. 10. 4X. ’RADICAL-’ . 

3  1P.D12.5.0P) 

END 


Function  RANGEO 


REAL *8  FUNCTION  RANGEO 
* - - - — 


LOGICAL*4  DEBUG 
INTEGER*4  I.  J.  K 

REAL*8  NU.  XP(3) .  R(3.2),  DX.  DY.  DZ.  DR 

REAL*8  Ml.  Al.  El.  INC1 .  NU1 .  ARGPAt .  LONAN1 .  XI (4) 

REAL*8  M2.  A2.  E2.  INC2,  NU2.  ARGPA2,  L0NAN2.  X2(4) 

REAL*8  HALFPI . PI . TWOPI . DEGRAD . RADDEG 

REAL*8  ER.  DU.  TU 

COMMON  /ADMIN/  DEBUG 

COMMON  /ASTRO/  ER.  DU.  TU 

COMMON  /COMPAR/  DX.  DY.  DZ.  DR 

COMMON  /CONST/  HALFPI.  PI.  TWOPI.  DEGRAD.  RADDEG 

COMMON  /SATl/  Ml.  Al .  El.  INCI .  NU1.  ARGPA1 .  LGNAN1 .  XI 

COMMON  /SAT2/  M2.  A2,  E2.  INC2.  NU2.  ARGPA2.  L0NAN2.  X2 


NU 1  -  NU(  Ml.  El) 


Compute  True  Anomaly  of  Sat  1 


Coaoute  the  Radius  of  Sat  1 


x:u:  *  a:  •  odd  -  E!*e:>  /  udo  *  ei 


Compute  the  position  of  Sat  1  in 
Perifocal  Coordinate  Frame 


DCCS (XU  1 ) ) 
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c 


c 

C 

C 

C 

C 

C 

C 


C 

c 

c 

c 


c 

c 

c 


* - * 

XPO)  -  X1(4)*DC0S(NU1) 

XP (2)  -  XI (4)*DSIN (NU1) 

XPO)  -  0D0 

* - * 

*  Compute  the  elements  of  the  * 

*  transformation  matrix  to  transform  * 

*  from  the  Perifocal  Coordinate  * 


*  Frame  to  the  Geocentric  Equatorial  * 

*  Frame.  * 

* — - — — - — - — - — * 

R(l.l)  -  DCQS ( LONAN 1 ) *DC0S (ARGPA1 ) - 
1  DSIN  (LONAN  1 )  *DSIN  (ARGPA1 )  *DC0S(INC1) 

RU.2>  »  -  DCQS  (LONAN  1 )  *DSIN  (ARGPA1 )  - 
1  DSIN  (LONAN  1 )  *DCQS  (ARGPA1 )  ♦DCQS  (INC1 ) 

R(2.!>  -  DSIN  (LONAN  1 )  *DCOS  (ARGPA1 )  ♦ 

1  DCQS (LONAN 1 ) *DSIN (ARGPA1) ♦DCOS ( INC1 ) 

R(2.2)  -  -  DSIN (LONAN l ) ♦DSIN (ARGPA1 ) * 

1  DCQS ( LONAN 1)*DC0S(ARGPA1)*DCQS( INC 1) 

R(3.i)  -  DSIN(ARGPA1)*DSIN(INC1) 

R(3.2)  -  DCOS(ARGPAl)  *DSIN  (INC1) 

IF  (DEBUG)  NRITE(*  .  1000)  (CR(I .  J) .  >1 ,2) .  I«1 .3) 

* — - - — - — 


*  Compute  the  position  of  Sat  1  in  * 

*  the  Geocentric  Equatorial  Frame.  * 

*  - — - — - — - — - * 

XI (1 >  «  R(1.1)*XP(1)  ♦  R(1.2)*XP(2) 

Xt (2)  •  R(2. l)eXP(l)  ♦  R(2.2)*XP(2) 

XI (3)  -  R(3. 1)*XP(1)  ♦  R(3.2)*XP(2) 


NU2  •  NU(  M2.  E2> 


Compute  the  True  Anomaly  of  Sat  2 


Comoute  the  Radius  of  Sat  2 


X2(4)  «  A2  *  (IDO  -  E2*E2)  /  (IDO  ♦  E2 


Compute  the  position  of  Sat  1  in 
Perifocal  Coordinate  Frame 


DCQS(NU2)) 


XP(1)  -  X2<4) *DCOS(NU2> 

XP(2)  «  X2(4)*DSIN(NU2) 

XPO)  «  ODC 

$ — - — - — - - - — - • 

*  Compute  the  elements  of  the  * 

♦  transformation  matrix  to  transform  ♦ 
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C  *  from  the  Perifocal  Coordinate  * 

C  *  Frame  to  the  Ceocentric  Equatorial  * 


C  *  Frame .  * 

C  * - - * 

R  ( 1 . 1 )  «  DCOS  ( L0NAN2)  ♦DCOS  (ARGPA2)  - 
1  DSIN  (L0NAN2)  ♦DSIN  (ARGPA2)  ♦DCOS(INC2) 

R(1 ,2)  -  -  DCOS (L0NAN2) ♦DSIN (ARGPA2)  - 
1  DSIN  CL0NAN2) ♦DCOS (ARGPA2) ♦DCOS ( INC2) 

R(2. 1)  -  DSIN (L0NAN2) ♦DCOS (ARGPA2)* 

1  DC0SCL0NAN2) ♦DSIN (ARGPA2) ♦DC0S(INC2) 

R(2,2)  -  -  DSIN (LONAN2) ♦DSIN (ARGPA2)* 

1  DCOS  (LONAN2)  ♦DCOS  (ARGPA2)  ♦DCOS  (INC2) 

R(3, 1)  -  DSIN (ARGPA2) ♦DSIN (INC2) 

R(3.2)  -  DCOS(ARGPA2)*DSIN CINC2) 

IF  (DEBUG)  WRITER , 1010)  (CRCI , J) . J-l .2) . 1*1 .3) 

C  ♦ . - . ♦ 

C  ♦  Compute  the  position  of  Sat  2  in  ♦ 

C  ♦  the  Geocentric  Equatorial  Frame.  ♦ 

C  ♦ - - ♦ 

X2C1)  -  R(1 . 1 ) ♦XP Cl)  ♦  R(1 .2) *XP(2) 

X2(2)  -  R(2. 1 ) ♦XPCl )  ♦  R(2.2)^XP(2) 

X2(3)  -  R(3. 1 ) ♦XP ( 1 )  ♦  R(3 .2) *XP(2) 


DX  ■  X2C1)  -  XI (1) 

DY  -  X2C2)  -  XI (2) 

DZ  •  X2(3)  -  XI (3) 

DR  -  DSQRT  CDX^DX  ♦  DY+DY  *DZ*DZ) 

IF  (DEBUG)  THEN 

WRITER . 1020)  NUl^RADDEG.  NU2^RADDEG 
WRITE(* . 1030)  (X2(I) . 1-1 .4) . (XI (I) . I«1 .4) . 

1  DX.  DY.  DZ.  DR 

ENDIF 

RANGE  -  DR 
RETURN 

1000  FORMAT  (///'  Sat  1  R-matrix  to  transform  from  Perifocal  \ 

1  'to  Geocentric  Equatorial  Frame.'/// 

2  (’  ' ,F15. 12.4X.F15. 12)//) 

1010  FOP‘”T  (///'  Sat  2  R-matrix  to  transform  from  Perifocal 

1  'to  Geocentric  Equatorial  Frame.'/// 

2  C  ' ,F15. 12.4X.F15. 12)//) 

1020  FORMAT  ('  Sat  1  True  Anomaly:  '.F16.10.4X. 

1  '  Sat  2  True  Anomaly:  ’.F16.10/) 

1030  FORMAT  (/T16. 'X ' .T34 . *Y 1 .T52. 'Z' .T66 . 'Magnitude '/ 

1  T9.' . 

2  * . 7 

3  '  Sat  2:  ' .F16.9.2X.F16.9.2X.F16.9.2X.F16.9// 

4  1  Sat  1:  ’.F16.9.2X.F16.9.2X.F16.9.2X.F16.9/ 

5  T9.' .  '. 
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6  ’ . ’/ 

7  ’  Delta:  ’ ,F16 .9 .2X.F16 .9 .2X.F16 .9 .2X.F16 .9//) 

END 

C  * - * 

C  *  * 

C  *  Function  NU  * 

C  *  ♦ 

C  * - - - — - * 

REAL*8  FUNCTION  NU C  M,  E  ) 

REAL*8  N.  E 

C  * - * 


LOG I CAL *4  DEBUG 
INTEGERS  I,  J.  K 

REAL*8  EANOM.  EANOM1 ,  EANDOT.  ERROR 
REAL*8  HALFPI . PI , TWOPI . DEGRAD . RADDEG 
REAL*8  ER.  DU.  TU 
COMMON  /ADMIN/  DEBUG 
COMMON  /ASTRO/  ER.  DU.  TU 

COMMON  /CONST/  HALFPI,  PI.  TWOPI .  DEGRAD.  RADDEG 

C  * - - - - — - * 

IF  <  E  .EQ.  0D0  )  THEN 
NU  -  M 

IF  (DEBUG)  WRITE  (*.1000)  NU 
RETURN 
ENDIF 
EANOM  -  M 

IF  (DEBUG)  WRITE  (*.1100)  M.  E.  EANOM 

I  -  0 

100  EANOM 1  -  EANOM 

1-1*1 

EANDOT  »  IDO  -  E  *  DCOS (EANOM 1) 

EANOM  -  EANOM  1  *  (M  -  EANOM  1  *  E*SIN(EANOMD) /EANDOT 
ERROR  -  DABS ( (EANOM 1  -  EANOM) /EANOM 1) 

IF  (DEBUG)  WRITE  (*.1200)  I.  EANOM 1 .  EANOM.  EANDOT.  ERROR 

IF  (ERROR  .GT.  5D-11)  GOTO  100 

NU  -  2D0  *  DATAN(  DSQRT((1D0  ♦  D/(1D0  -  D)  * 

1  DTAN(5D-1  *  EANOM)) 

IF  (NU.LT.O)  NU  -  TWOPI  ♦  NU 
RETURN 

1000  FORMAT  (/’  This  is  a  circular  orbit,  so  NU  -  M.  NU»\ 

1  F12.10) 

1100  FORMAT  (//T4, ‘Mean’. T35. ’Eccentric’/ 

1  T3. ’Anomaly ’ .T18. ’Eccentricity ’ .T36. ’Anomally'/ 

2  T2.  ’ . - . 7 

3  1X.3(F12. 10.  4X)/// 

4  T18. ’Old’ .T35. ’New’/T16. ’Eccentric ’ .T33. ’Eccentric’/ 

5  ’  Iteration* .T17. ’Anomaly ’ .T34. ’Anomaly* .T50. ’dM/dE’ . 
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I 

^  6  T66, ’Error’/ ’  . 

17  ’ . ’) 

1200  FORMAT  (2X.  17,  5X,  F12.10,  4X,  F12.10,  4X,  IP.  D12.5.  4X. 
1  D12.5,  0P> 

|  END 

I 

I 

I 

I 

E 

I 

I 


♦include  <stdio.h> 


/* - */ 

/*  */ 

/*  rndini  -  this  routine  initializes  the  random  number  */ 
/*  generator.  */ 

/*  */ 

/* - */ 


rndini 0 

< 

srand(l) ; 


> 

/* - ♦/ 

/*  */ 

/*  random  -  this  routine  returns  a  double  precision  */ 

/*  uniform  distributed  random  variable  between  */ 

/*  0  and  1 .  */ 

/*  ♦/ 

/* - - */ 


double  *ranaom(> 

int  rand  O  ; 
double  r; 
r  «  randO ; 
r  *«  32768.0; 
r  randO  ; 
r  /-  32768.0  *  32768.0; 
return  <&r) ; 


00 
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CGPCA  -  Circular  Orbit  Probability  of  Close  Approach 


PROGRAM  COPCA 
CHARACTER  YESNO*4 
LOGICAL*4  DEBUG 
INTEGER*4  ONE,  TWO 

REAL*8  Ml .  Ai,  El .  INC! t  NU1 .  ARGPA1 ,  LONAN1 .  XI (4) 

REAL*8  M2.  A2,  E2.  INC2.  NU2,  ARGPA2.  L0NAN2.  X2C4) 

REAL*8  ANGDTH .  APCA,  A2PCA,  APXPCA,  CADTH 

RE4L*8  DTH,  INCR 

REAL*8  M.  OPSEP.  PCA.  PDCA 

REAL*8  PDCAM1 .  PDCAM2.  PDCAM3.  SADTH.  STHETA.  THETA 

REAL*8  HALFPI . PI . TWOPI . DEGRAD . RADDEG 

REAL*8  ER.  DU.  TU 

COMMON  /ADMIN/  DEBUG 

COMMON  /ASTRO/  ER.  DU.  TU 

COMMON  /CONST/  HALFPI.  PI.  TWOPI,  DEGRAD.  RADDEG 
COMMON  /SAT1/  Ml.  A1 .  El.  INC1 .  NU1.  ARGPA1 ,  LONAN1 .  XI 

COMMON  /SAT2/  M2.  A2.  E2.  INC2.  NU2.  ARGPA2.  L0NAN2.  X2 

COMMON  /TRANS/  SADTH.  CADTH.  STHETA 
PARAMETER  CONE-1.  TWO-2) 

ER  «  6378. 145D0 
DU  -  6378. 145D0 
TU  -  806.81 18744D0 
PI-DACOS(-iDO) 

HALFPI-PI/2D0 

TW0PI-2D0*PI 

DEGRAD-PI/180D0 

RADDEG-180D0/PI 


C  * - * 

C  *  Incut  Orbital  Elements  of  Sat  1  * 

C  * . * 


WRITE (♦.1000)  ONE 
READ  (♦.1010)  A1 
WRITEC* , 1020)  ONE 
READ  (* .1010)  El 
WRITE (♦.1030)  ONE 
READ  (♦.1010)  INC1 
WRITER.  1040)  ONE 
READ  (♦.1010)  ARGPA1 
WRITEC*. 1050)  ONE 
READ  (*.1010)  LON AN! 

C  * . * 

C  ♦  Input  Orbital  Elements  of  Sat  2  ♦ 
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c  * - — — * 

WRITE (*. 1000)  TWO 
READ  (*,1010)  A2 
WRITE (*, 1020)  TWO 
READ  (*,1010)  E2 
WRITE(* . 1030)  TWO 
READ  (*,1010)  INC2 
WRITE(*, 1040)  TWO 
READ  (*.1010)  ARGPA2 
WRITE (*, 1050)  TWO 
READ  (*.1010)  L0NAN2 


C  * — - - — - * 

C  *  Input  Simulation  Limits  * 

C  * - * 


WRITE (*.1060) 

READ  (*,1010)  DTK 
WRITE (*.1070) 

READ  (*.1080)  ITER 
I NCR  -  TWOPI  /  DBLE(ITER) 

WRITE(*, 1090) 

READ  (*.1100)  YESNO 

IF  (YESNOO  :  1)  .EQ.  ’Y  *  .OR.  YESNOO : 1) .EQ. *y’)  THEN 
DEBUG  -  .TRUE. 

ELSE 

DEBUG  -  .FALSE. 

ENDIF 

C  * - * 

C  *  Convert  the  input  angles  from  deg  * 

C  *  to  rad  * 

C  * - ♦ 

INC1  «  INC1  *  DEGRAD 
ARGPA1  -  ARGPA1  *  DEGRAD 
LONAN1  -  LONAN1  *  DEGRAD 
INC2  -  INC2  *  DEGRAD 
ARGPA2  -  ARGPA2  *  DEGRAD 
L0NAN2  -  L0NAN2  *  DEGRAD 

C  * - - ♦ 

C  *  Initialize  Probability  of  Close  • 

C  ♦  Approach  * 

C  * - * 

PCA  -  0D0 

C  * - - - * 

C  *  Compute  Angular  Separation  Between  * 

C  *  The  Orbital  Planes  of  Sat  1  and  * 

C  *  Sat  2  * 

C  * - - - * 

THETA  -  OPSEPO 
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STHETA  -  DSIN (THETA) 

C  * - * 

C  *  Conpute  Angular  Distance  Threshold  * 

C  * - - * 

ANGDTH  -  (A1*A1  ♦  A2*A2  -  DTH*DTH)/(2D0*A1*A2) 
IF  (DABS (ANGDTH)  .GT.  IDO)  THEN 
IF  (ANGDTH  .LT.  -IDO)  THEN 
ANGDTH  -  PI 
ELSE 

ANGDTH  -  ODO 
END  IF 
ELSE 

ANGDTH  -  DACOS (ANGDTH) 

ENDIF 

SADTH  -  DSIN (ANGDTH) 

GADTH  -  DCQS (ANGDTH) 


C  ♦ - - - ♦ 

C  *  Compute  Approximate  Probability  of  * 

C  *  Close  Approach  (Using  Both  Methods)* 

C  ♦ - * 

APCA  -  APXPCA (ANGDTH) 

C  ♦ - * 

C  *  Start  Numerical  Integration  * 

C  * - * 


IF  (ANGDTH  .GE.  THETA)  THEN 
M2  -  HALFPI 
ELSE 

M2  •  DASIN(SADTH  /  STHETA) 

ENDIF 
Ml  »  -M2 
IF  (DEBUG)  THEN 
WRITER.  11 10)  Ml.  M2 
M2  •  Ml  ♦  IDI  *  I NCR 
ENDIF 
M  ■  Ml 

C  * . . — * 

C  *  Adjust  step  size  to  make  the  * 

C  *  Integration  range  an  Integer  * 

C  *  number  of  step  size’s.  * 

C  — — — - - * 

IF  (  M2  ,GT.  Ml)  I NCR  ■  (M2-M1)/DBLE(IDNINT((M2-M1)/INCR)> 
M2  -  M2  -  INCR/2D0 
PDCAM1  -  PDCA(M) 

100  PDCAM2  •  PDCACM+ .5D0*INCR) 

PDCAM3  »  PDCA(M*INCR) 

PCA  ■  PCA  ♦  INCR/6D0  *  (PDCAM1  ♦  4D0*PDCAM2  ♦  PDCAM3) 

IF  (DEBUG)  NRITE(* . 1 120)  M.  PDCAM1 .  PDCAM2.  PDCAM3.  PCA 
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PDCAM1  «  PDCAM3 
M  ■  M  ♦  I  NCR 
IF  (N  .LE.  N2)  GOTO  100 
PCA  -  2D0  *  PCA 

WRITE(* , 1 130)  Al.  A2.  FI.  E2.  INC1*RADDEG.  INC2*RADDEG, 

1  ARGPA1*RADDEG,  ARGPA2*RADDEG ,  LONAN 1 *RADDEG , 

2  L0NAN2*RADDEG,  INCR*1D6.  DTH,  ANGDTH*RADDEG , 

3  THETA*RADDEG .  APCA,  PCA 

C  ♦ - ♦ 

1000  FORMAT  (/'  Input  the  semi-major  axis  of  Sat  *.I1, 

1  '  in  km:  ’\) 

1010  FORMAT  (F16.12) 

1020  FORMAT  ('  Input  the  eccentricity  of  Sat  Ml.*:  *\) 

1030  FORMAT  <*  Input  the  inclinition  of  Sat  Ml, 

1  *  in  degrees:  *\) 

1040  FORMAT  (*  Input  the  argument  of  perapsis  of  Sat  '.II. 

1  *  in  degrees:  *\) 

1050  FORMAT  C*  Input  the  longitude  of  the  ascending  node  *. 

1  'of  Sat  Ml.'  in  degrees:  ’\) 

1060  FORMAT  (/'  Input  distance  threshold  in  km:  *\) 

1070  FORMAT  ('  Input  the  desired  number  of  iterations  per  radian  *. 

1  '<7  digits  max) :  '\> 

1080  FORMAT  (17) 

1090  FORMAT  (’  Run  COPCA  in  DEBUG  mode  (Y/N)?  ’\) 

1100  FORMAT  (A4) 

1110  FORMAT  <//'  Ml»* . F13. 10.4X, 'M2*  * ,F13. 10/) 

1120  FORMAT  ('  M-'.F16.10/'  PDCAM1-’ .F16. 10/'  PDCAM2-* .F16 ,10/ 

1  '  PDCAM3*' ,F16. 10/'  PCA  «\F16.10///) 

1130  FORMAT  C/T21 .  'Sat  l'.T41.'Sat  2VT16.' . *. 

1  T36. ' . '//'  a  (km) '  .T16.F16. 10. 4X. 

2  F16.10//'  e* .T16.F16. 10.4X.F16. 10// 

3  *  inc  (deg) ’ .T16.F16. 10.4X.F16. 10// 

4  '  argument  of ’ .T16.F16. 10.4X.F16. 10/'  perigee  (deg)'// 

5  '  long  of  asc' .T16.F16. 10.4X.F16. 10/'  node  (deg)'/// 

6  ’  Mean  Anomaly  Iteration  Step  Size  :  '. 

7  F16.10,’  micro-radians’// 

8  '  Distance  Threshold  for  Close  Approach: 

9  F16. 10. '  km’// 

A  *  Angular  Threshold  for  Close  Approach  :  ’. 

B  F16.10.'  deg'// 

C  '  Angular  Separation  of  Orbital  Planes  :  '. 

D  F16. 10. *  deg'// 

E  *  Approximate  Probability  of  Close  Approach:  *. 

F  F9  6// 

G  '  Computed  Probability  of  Close  Approach  :  ’.F9.6//) 

END 

C  * - ♦ 
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Function  OPSEPO 


* 


♦ 


♦ 

♦ 


* 

♦ 


♦ - - - * 

REAL*8  FUNCTION  OPSEP 

* — — - - — * 

LOGICAL^  DEBUG 
REAL*8  DLONAN,  THETA 

REAL*8  Ml.  Al.  El.  INC1 .  NU1 .  ARGPA1 .  L0NAN1 .  XI (4) 

REAL*8  M2.  A2.  E2.  INC2.  NU2.  ARGPA2,  LONAN2.  X2C4) 

REAL*8  HALFPI . PI . TWOPI . DEGRAD . RADDEG 
REAL*8  ER.  DU.  TU 
COMMON  /ADMIN/  DEBUG 
COMMON  /ASTRO/  ER.  DU.  TU 

COMMON  /CONST/  HALFPI.  PI.  TWOPI.  DEGRAD.  RADDEG 
COMMON  /SAT1/  Ml.  Al .  El.  INC1 .  NU1 .  ARGPA1 .  LONAN1.  XI 

COMMON  /SAT2/  M2.  A2.  E2.  INC2.  NU2.  ARGPA2.  L0NAN2.  X2 


DLONAN  -  L0NAN2  -  LONAN1 
IF  (DABS (DLONAN)  .GT.  PI) 

1  DLONAN  -  DLONAN  -  DSIGN (TWOPI  .DLONAN) 
THETA  -  DAC0S(DC0S(INC1)*DC0S(INC2) 

1  ♦  DSIN(INC1)*D5IN(INC2)*DC0S(DL0NAN)) 

IF  (THETA  .GT.  HALFPI)  THETA  -  PI  -  THETA 
OPSEP  -  THETA 
RETURN 
END 


Function  PDCAO 


REAL*8  FUNCTION  PDCA(M) 

C  * . . * 

LOCI CAL *4  DEBUG 

REAL*8  ARC.  M.  SM.  STHETA,  SADTH.  CADTH 

REAL+8  Ml.  Al.  El.  INC1 .  NU1.  ARGPA1 .  LON AN 1 .  XI (4) 

REAL*8  M2.  A2.  E2.  INC2.  NU2,  ARGPA2.  LONAN2.  X2(4) 

REAL*8  HALFPI . PI . TWOPI . DEGRAD . RADDEG 
REAL*8  ER.  DU.  TU 
COMMON  /ADMIN/  DEBUG 
COMMON  /ASTRO/  ER.  DU.  TU 

COMMON  /CONST/  HALFPI.  PI.  TWOPI.  DEGRAD.  RADDEG 
COMMON  /SATl/  Ml.  Al .  El.  INC!.  NU1 .  ARCPA1 .  LON AN! .  XI 

COMMON  /SAT2/  M2.  A2.  E2.  INC2.  NU2.  ARGPA2.  L0NAN2.  X2 

COMMON  /TRANS/  SADTH.  CADTH.  STHETA 
C  * - ♦ 
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SN  -  DSIMC  M  ) 

ARG  «  IDO  -  STHETA*STHETA  ♦  SM*SN 
IF  (ARG  .EQ.  ODO)  THEN 
IF  (CADTH  .GE.  ODO)  THEN 
PDCA  -  ODO 
RETURN 

&LDJL 

PDCA  ■  5D-1  /  ,'I 
RETURN 
ENDIF 
END  IF 

ARG  •  CADTH  /  DSQRT(ARG) 

IF  (DABS (ARC)  .GT.  IDO)  THEN 
IF  (ARG  .GT.  IDO)  THEN 
PDCA  ■  ODO 
RETURN 
ELSE 

PDCA  ■  5D-1  /  PI 
RETURN 
ENDIF 
ENDIF 

PDCA  ■  DACOS (ARG) / (2D0  *  PI  *  PI) 

RETURN 

END 


Function  APXPCAO 


C 


C 


REAL*8  FUNCTION  APXPCA(ANGDTH) 
REALMS  ANGDTH 


• - • 

LOCICAL'4  DEBUG 


REAL'S  STHETA.  SADTH.  CADTH 

REAL'8  Ni .  Ai.  El.  INC1 .  NU1.  ARCPAi.  LON AN 1 .  XI (4) 
REAL'S  N2.  A2.  E2.  INC2.  NU2.  ARGPA2.  LONAN2.  X2U) 
REAL'8  HALFPI .PI .TNOPI . DEGRAD. RADDEC 


REAL'S  ER.  DU.  TU 
COMMON  /ADMIN/  DEBUG 


COMMON  /ASTRO/  ER.  DU.  TU 

COMMON  /CONST/  HALFPI.  PI.  TNOPI.  DEGRAD  RADDEC 


COMMON  /SMI/ 
CO IOON  /SAT2/ 


Ml.  Al.  El.  INC! .  NUl .  ARCPAI.  LON AN 1 .  XI 
M2.  A2.  E2.  INC2.  NU2.  ARCPA2.  L0NAN2.  X2 


COMMON  /TRANS/  SADTH.  CADTH.  STHETA 
•--- - ..... - • 


IF  (STHETA  .LT.  ID-1)  THEN 
APXPCA  •  -IDO 
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ELSE 

APXPCA  -  AMGDTH  *  AXGDTH  /  (2D0  *  PI  ♦  STHETA) 
IF  (APXPCA  .GT.  IDO)  APXPCA  »  -IDO 
ENDIF 
RETURN 
END 


o  o  o  o  o 


EOPCA  -  Elliptical  Orbit  Probability  of  Close  Approach 


PROGRAM  EOPCA 
CHARACTER  TESH04 


LOGICAL*4  DEBUG 
IITECER*4  ONE.  TNO 
REAL+8  N 

REAL*8  BOUND  1 .  B0UND2.  DTH.  DSQR,  PDTH.  PR.  NUU.  NU12. 
1  NU21 .  NU22 

REAL*8  Ml .  A1 .  El .  INC1 .  NU1 .  ARGPA1 .  LON AN 1 .  Xl <45 
REAL*8  M2.  A2.  E2.  INC2.  NU2.  ARGPA2.  LONAN2.  X2(4> 
REAL*8  SINC1 .  CINC1 .  SNU1 .  CNU1 .  SARGP1 .  CARGPl . 

1  StXANl .  CLNAN1 .  BETA! 

REAL’S  SINC2.  CINC2.  SNU2.  CNU2.  SARGP2.  CARGP2. 

1  SLNAN2.  CLNAN2.  BETA2 

REAL’S  Rll .  R12.  R13.  R21 .  R22.  R23.  R31.  R32.  R33 

REAL’S  SU.  S12.  S13.  S21.  S22.  S23.  S31.  S32.  S33 

REAL’S  AFCA 

REAL’8  INCR.  LIMIT.  NU 

REAL*8  PCA.  PDCA 


REAL’S  PDCAN1 .  PDCAM2.  PDCAM3 
REAL’S  HALFPI . PI . TNOPI . DEGRAD . RADDEG 


REAL’S  ER.  DU.  TU 
IN  /ADMIN/  DEBUG 


CONNOR  /ASTRO/  ER.  DU.  TU 

CONNOR  /BOURD/  BOURD 1 .  BOURD2.  DTH.  DSQR.  PDTH.  PR.  MU  11. 

1  MU12.  MU21 .  MU22 

CONNOR  /CORST/  HALFPI.  PI.  TROPI.  DEGRAD.  RADDEG 
CONNOR  /SATI/  Ml.  A1 .  El.  IRC1 .  RUl.  ARGPA1 .  LORAR1.  XI. 
t  SIRC1.  CIRC1.  SRU1.  CRU1.  SARGP1 .  CARGPl. 

2  SLRAR1 .  CLRAR1.  BETA! 

CONNOR  /SAT2/  N2.  A2.  E2.  IRC2.  RU2.  ARGPA2.  LORAR2.  X2. 

1  SIRC2.  CIRC2.  SRU2.  CRU2.  SARGP2.  CARCP2. 

2  SLRAR2.  CLRAR2.  BETA2 


CONNOR  /TRAMS/  Rll.  R12.  R13.  R21 .  R22.  R23.  R31 .  R32. 
1  R33.  Sll.  $12.  S13.  S21 .  S22.  S23.  S31 .  S32.  S33 


PARAMETER  (ORE-1 .  TRO-2) 
ER  •  6373. 145DO 


DU  •  6378.14SD0 
TU  •  806.81 18744DO 


PI  •  daxsc-idc: 
HALFPI  •  PI/2D0 


TROPI  •  2D0*PI 
DEGRAD  •  PI/180D0 
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RADDEG  «  180D0/PI 

* - * 

*  Input  Orbital  Elements  of  Sat  1  * 

♦  - *  • 


WRITE (*,1000)  ONE 
READ  (*,1010)  A1 
WRITE (*,1020)  ONE 
READ  (*.1010)  El 
WRITE (*,1030)  ONE 
READ  (*.1010)  INC1 
WRITE (*,1040)  ONE 
READ  (*,1010)  ARGPA1 
WRITE(* , 1050)  ONE 
READ  (*.1010)  LONAN1 


♦ - * 

*  Input  Orbital  Elements  of  Sat  2  * 

*  - * - - - * 


WRITE(* , 1000)  TWO 
READ  (*.1010)  A2 
WRITE (* , 1020)  TWO 
READ  (*.1010)  E2 
WRITE(*, 1030)  TWO 
READ  (*,1010)  INC2 
WRITE(*. 1040)  TWO 
READ  (*,1010)  ARGPA2 
WRITE (*,1050)  TWO 
READ  (*.1010)  L0NAN2 


* - — - - * 

*  Input  Simulation  Limits  * 

*  - - - - — — - - -* 


WRITE(* . 1060) 

READ  (*.1010)  DTH 
WRITE (*,1070) 

READ  (*.1080)  ITER 
I NCR  «  TWOPI  /  DBLE(ITER) 

WRITE(* . 1090) 

READ  (*.1100)  YESNO 

IF  (YESNO(l : 1) .EQ. *Y* .OR. YESNO(l : 1) .EQ. *y *>  THEN 
DEBUG  «  .TRUE. 


ELSE 

DEBUG  -  .FALSE. 

END  IF 

* - - * 

*  Convert  the  input  angles  from  deg  * 

*  to  rad  * 

*  - - * 


INC1  •  INC1  *  DEGRAD 
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ARGPA1  *  ARGPA1  *  DEGRAD 
L0NAN1  -  LONAN1  *  DEGRAD 
INC2  -  INC2  *  DEGRAD 
ARGPA2  *  ARGPA2  *  DEGRAD 
L0NAN2  *  L0NAN2  *  DEGRAD 

C  * . * 

C  *  Compute  the  elements  of  the  * 

C  *  transformation  matrix  to  transform  * 

C  *  from  the  Perifocal  Coordinate  * 

C  *  Frame  to  the  Geocentric  Equatorial  * 


C  *  Frame .  * 

C  * . * 


SARGP1  -  DSINCARGPA1) 
CARGP1  -  DCOS(ARGPAl) 
SINC1  -  DSIN(INCl) 
CINC1  -  DCOS(INCl) 
SLNAN 1  -  DSIN(LONANl) 
CLNAN 1  *  DCOS(LONANl) 


BETA1  -  DSQRT(<1D0-E1)/(1D0+E1)) 

C  * - * 


Rll  -  CLNAN 1 +CARGP 1 -SLNAN 1 *SARGP 1 *CI NC 1 
R12  -  -CLNAN 1 *SARGP 1 -SLNAN 1 ♦CARGP 1 *CI NC 1 
R13  -  0D0 

R21  -  SLNAN 1 ♦CARGP1 +CLNAN 1 *SARGP1 *CINC1 

R22  -  -SLNAN 1 ♦SARGP1 +CLNAN 1 +CARGP1 ♦CINCl 

R23  -  ODO 

R31  -  SARGP1*SINC1 

R32  -  CARGP1*SINC1 

R33  -  ODO 

IF  (DEBUG)  NRITE<* .1110)  BETA1 ,R11 ,R12,R21 ,R22,R3i ,R32 

Q  * - - * 

C  *  Compute  the  elements  of  the  * 

C  *  transformation  matrix  to  transform  * 

C  *  from  the  Perifocal  Coordinate  * 

C  *  Frame  to  the  Geocentric  Equatorial  * 

C  ♦  Frame .  * 

C  * - — * 

SARGP2  -  DSINCARGPA2) 

CARGP2  *  DCOS (ARGPA2) 

SINC2  -  DSINCINC2) 

CINC2  «  DC0SCINC2) 

SLNAN2  -  DSINCL0NAN2) 

CLNAN2  *  DC0SCL0NAN2) 

BETA2  -  D5QRT((1D0-E2)/(1D0*E2>) 

C  * . . * 

511  -  CLNAN2*CARGP2-SLNAN2*SARGP2*CINC2 

512  -  SLNAN2*CARGP2*CLNAN2*SARGP2*CINC2 
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S13  ■  SARGP2*SINC2 

521  -  -CLNAN2+SARGP2-SLNAN2*CARGP2*CI NC2 

522  »  -SLNAN2*SARGP2+CLNAN2*CARGP2*CINC2 

523  »  CARGP2*SINC2 

531  -  SLNAN2*SINC2 

532  «  -CLNAN2*SINC2 

533  «  CINC2 

IF  (DEBUG)  WRITE (*, 1120)  BETA2.S11,  S12,  S13,  S21 .  S22. 


1  S23.  S31,  S32,  S33 

C  * - * 

C  *  Initialize  Probability  of  Close  * 

C  *  Approach  * 

C  * - - - . * 


PCA  »  0D0 

C  * - * 

C  *  Start  Numerical  Integration  * 

C  * - - * 

NU  -  -PI 

LIMIT  -  PI  -  INCR/2D0 
PDCAM1  -  PDCA(NU) 

100  PDCAM2  -  PDCA(NU+.5D0*INCR) 

PDCAM3  -  PDCACNU+INCR) 


PCA  -  PCA  ♦  INCR/6D0  *  CPDCAM1  ♦  4D0*PDCAM2  ♦  PDCAN3) 

NU  -  NU  ♦  INCR 

IF  (DEBUG)  WRITE (* ,1130)  NU,  PDCAM1 .  PDCAM2.  PDCAN3.  PCA 
PDCAM1  -  PDCAM3 
DEBUG  -  .FALSE. 

IF  (NU  .LT.  LIMIT)  GOTO  100 

WRITE(* , 1140)  Al.  A2,  El.  E2.  INC1*RADDEG,  INC2*RADDEG, 

1  ARGPA1*RADDEG.  ARGPA2*RADDEC .  LONAN 1 *RADDEG . 

2  L0NAN2*RADDEG .  INCRMD6.  DTK.  PCA 


C  * - 

1000  FORMAT 
1 

1010  FORMAT 
1020  FORMAT 
1030  FORMAT 
1 

1040  FORMAT 
1 

1050  FORMAT 
1 

1060  FORMAT 
1070  FORMAT 
1 

1080  FORMAT 
1090  FORMAT 


(/'  Input  the  semi-major  axis  of  Sat  Ml, 

1  in  km:  *\) 

(F16.12) 

('  Input  the  eccentricity  of  Sat  Ml.*:  '\) 

('  Input  the  inclinition  of  Sat  Ml, 

*  in  degrees:  '\) 

(’  Input  the  argument  of  perapsis  of  Sat  Ml. 
*  in  degrees:  ’\) 

(’  Input  the  longitude  of  the  ascending  node 
•of  Sat  Ml.'  in  degrees:  '\) 

(/'  Input  distance  threshold  in  km:  *\) 

(’  Input  the  desired  number  of  iterations  per 
'radian  (7  digits  max) :  '\) 

(17) 

('  Run  EOPCA  in  DEBUG  mode  (Y/N)?  *\) 


♦ 
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1100  FORMAT  CA4) 

1110  FORMAT  (///*  Betal  -  ’.F15.12// 

1  *  Sat  1  R-matrix  to  transform  from  Perifocal  to  ’ , 

2  *  Geocentric  Equatorial  Frame.*// 

3  (*  * ,F15. 12.4X.F15. 12.4X. *  0.0*)/) 

1120  FORMAT  (//*  Beta2  *  ’.F15.12// 

1  *  Sat  2  S-matrix  to  transform  from  Geocentric  *, 

2  ’Equatorial  Frame  to  Sat  2**s  Perifocal  Frame.*// 

3  (*  * ,F15. 12.4X.F15. 12.4X.F15. 12)/) 

1130  FORMAT  (*  NU  -’.F16.10/’  PDCAM1-* ,F16. 10/ *  PDCAM2* * . 

1  F16.10/*  PDCAM3*’ . FI 6 . 1 0/ *  PCA  «’ . F16 . 10///) 

1140  FORMAT  C/T21 ,  *Sat  l’,T41.’Sat  2*/Tl6.* - *. 

1  T36 .  * - *//*  a  (km)  *  .T16.F16. 10.4X, 

2  F16.10//*  e* .T16.F16. 10.4X.F16 . 10// 

3  *  inc  (deg) * .T16.F16. 10.4X.F16. 10// 

4  *  argument  of  * .T16.F16.10.4X.F16.1C/*  perigee  (deg)*// 

5  *  long  of  asc *.T16.F16.10,4X,F16.10/*  node  (deg)*/// 

6  *  True  Anomaly  Iteration  Step  Size  :  *, 

7  F16.10,’  micro-radians*// 

8  ’  Distance  Threshold  for  Close  Approach: 

9  F16.10.’  km’// 

A  *  Computed  Probability  of  Close  Approach  :  *,F9.6//) 


END 

C  * - * 

C  *  * 

C  *  Subroutine  ORBBND  * 

C  *  * 

C  * . * 

SUBROUTINE  ORBBND () 

C  * — - — - — — * 

LOG  I  CALM  DEBUG.  INSIDE 
INTEGERS  I.  N.  NORDER 


REAL*8  A.  B.  C.  D.  E.  F.  G.  H.  K.  SAT2P 
REAL*8  PI.  P2.  P3.  P4,  P5.  P(4) 

REAL*8  SOLUTION (4).  RANK (4) 

REAL*8  ANGLE.  DIST,  TEMPI.  TEMP2 
REAL*8  NU.  X.  Y.  Z 

REAL*8  BOUND 1 .  B0UND2.  DTH.  DSQR.  PDTH,  PR.  NU11.  NU12. 
1  NU21 .  NU22 

REAL*8  Ml .  A1 .  El .  INC1 .  NU1 .  ARGPA1 .  LONAN1 .  XI (4) . 

1  SINC1.  CINC1.  SNU1.  CNU1 .  SARGP1 .  CARGP1 . 

2  SLNAN1 .  CLNAN1.  BETAl 

REAL*8  M2.  A2.  E2.  INC2,  NU2.  ARGPA2.  L0NAN2.  X2(4) . 

1  SINC2.  CINC2.  SNU2.  CNU2.  SARGP2 .  CARGP2 . 

2  SLNAN2.  CLNAN2.  BETA2 

REAL*8  Rll.  R12.  R13,  R21 .  R22.  R23.  R31 .  R32.  R33. 

1  Sll,  S12.  S13.  S21.  S22.  S23.  S31 .  S32.  S33 
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REAL*8  HALFPI , PI , TWOPI , DEGRAD , RADDEG 
REAL*8  ER,  DU,  TU 
COMPLEX* 16  ANSWER (4) 

COMMON  /ADMIN/  DEBUG 
COMMON  /ASTRO/  ER,  DU,  TU 

COMMON  /BOUND/  BOUND 1 ,  B0UND2 ,  DTH,  DSQR,  PDTH,  PR,  NU11, 
1  NU12,  NU21 ,  NU22 

COMMON  /CONST/  HALFPI.  PI,  TWOPI,  DEGRAD,  RADDEG 
COMMON  /SAT1/  Ml,  A1 ,  El,  INC1 ,  NU1 ,  ARGPA1 ,  LONAN1,  XI, 

1  SINC1 ,  CINC1 ,  SNU1,  CNU1 ,  SARGP1 ,  CARGP1 , 

2  SLNAN1,  CLNAN1 ,  BETA1 

COMMON  /SAT2/  M2,  A2,  E2,  INC2,  NU2,  ARGPA2,  L0NAN2 ,  X2, 

1  SINC2,  CINC2,  SNU2,  CNU2,  SARGP2,  CARGP2 , 

2  SLNAN2 ,  CLNAN2 .  BETA2 

COMMON  /TRANS/  Rll,  R12,  R13,  R21 ,  R22,  R23,  R31 ,  R32, 

1  R33 ,  Sll,  S12,  S13 ,  S21,  S22,  S23,  S31,  S32,  S33 


* - * 

*  Initialize  Arrays  to  zero  * 

* — — - - — - — - * 


DO  100  1-1,4 
ANSWER(I)-(0D0 , 0D0) 

P(I)-0D0 
RANKCD-ODO 
SOLUTION  (I)-ODO 
100  CONTINUE 

C  * - * 

C  *  Find  the  terms  needed  to  compute  * 

C  *  the  polynomial  coefficients  of  the  * 

C  *  polynomial  used  to  find  NU11  thru  * 

C  *  NU22.  Use  2nd  order  polynomial  * 

C  *  solution  if  XI (l)  or  XI (2)  are  * 

C  *  very  small.  * 

C  *— - — - - - - — — * 

SAT2P  -  A2*(1D0  -  E?*E2) 

A  •  E2  *  ((DSQR  -  PR* PR) *E2  *  2D0*X1 (1)*SAT2P) 

B  -  2D0* ( (DSQR  -  PR*PR)*E2  ♦  X1(1)*SAT2P> 

C  -  DSQR  -  PR* PR  -  SAT2P*SAT2P 

K  -  4D0*X1 (2) *X1 (2) *SAT2P*SAT2P 

IF  (DEBUG)  WRITE(* ,900)  DSQR.  PR,  E2.  SAT2P,  XI (1),  A,  B, 


1  C,  K 

C  * . * 


IF  (DABS (XI (2)) .GT. ID-8  .AND.  PR , NE . DABS (X 1 (1)))  THEN 
PI  -  A*A  *  K*E2*E2 
P2  -  2D0* (A*B  ♦  K*E2) 

P3  -  B*B  *  2D0*A*C  *  K*(1D0-E2*E2) 

P4  «  2D0*(B*C  -  K*E2) 
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P5  *  C*C  -  K 

* - * 

*  Divide  the  polynomial  by  the  * 

*  highest  order  coefficient  * 

*  - * 

NORDER  «  0 
IF  CP1.NE.0D0)  THEN 
♦ - * 

*  The  polynomial  is  4th  order  * 

*  - * 


NORDER  »  4 
PCD  «  P2/P1 
P(2)  «  P3/P1 
P(3)  «  P4/P1 
P (4)  «  P5/P1 

CALL  QUART ICC  P,  ANSWER) 

ELSE  IF  (P2.NE.0D0)  THEN 

♦ - - - * 

*  The  polynomial  is  3rd  order  * 

*  - - - -* 

NORDER  -  3 

PCI)  *  P3/P2 
PC2)  -  P4/P2 
PC3)  -  P5/P2 
CALL  CUBIC C  P.  ANSWER) 

ELSE  IF  (P3.NE.0D0)  THEN 

* - — — - - - - - * 

♦  The  polynomial  is  2nd  order  * 

*  - * 

NORDER  -  2 
PCI)  -  P4/P3 
PC2)  «  P5/P3 

CALL  QUADRATIC C  P,  ANSWER) 

ENDIF 

ELSE 

IF  (A.NE.ODO)  THEN 
NORDER  -  4 
PCD  -  B/A 
PC2)  -  C/A 

CALL  QUADRATIC C  P,  ANSWER) 

ANSWER  C3)  -  ANSWER  CD 
ANSWER C4)  -  ANSWER C2) 

ELSE 

NORDER  «  2 
ANSWER  CD  -  -C  /  B 
ANSWER  C2)  «  ANSWER  CD 
ENDIF 
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ENDIF 

IF  (DEBUG)  THEN 

WRITE(*,1000)  NORDER.  PI,  P2,  P3,  P4,  P5 
DO  200  1-1,4 

WRITE (* ,1010)  I,  P(I).  I,  DREAL (ANSWER (I) ) , 

1  DIMAG  (ANSWER(I)) 

200  CONTINUE 
ENDIF 
N  -  0 

DO  300  I-l.NORDER 

IF  (DIMAG (ANSWER ( I) ).EQ.0D0)  THEN 

IF  ((DABS (DREAL (ANSWER(I))) -IDO) .LE.1D-10)  THEN 
N  -  N  ♦  1 

SOLUTION (N)  «  DREAL (ANSWER (I)) 

IF  (DABS (SOLUTION (N)).GT. IDO)  THEN 
IF  (SOLUTION (N) .GT.ODO)  THEN 
SOLUTION (N)  -  0D0 
ELSE 

SOLUTION (N)  -  PI 
ENDIF 
ELSE 

SOLUTION (N)  -  DACOS (SOLUTION (N)) 

ENDIF 
ENDIF 
ENDIF 
300  CONTINUE 

IF  (DEBUG)  WRITER.  1020)  N,  (I  ,RADDEG*SOLUTION(I) .  1-1  ,N) 
IF  (N.GE.2)  THEN 
RANK(l)  -  SOLUTION (1) 

DO  500  1-2. N 
TEMP2  »  SOLUTION (I) 

DO  400  >1.1-1 
IF  (TEMP2.LT.RANK(J))  THEN 
TEMPI  -  RANK(J) 

RANKuf)  -  TEMP2 
TEMP2  •  TEMPI 
ENDIF 

400  CONTINUE 

RANK (I)  -  TEMP2 
500  CONTINUE 
ENDIF 

IF  (DEBUG)  WRITER.  1030)  N.  (I  ,RADDEG*RANK(I)  .1-1  .N) 


C  ♦ - - * 

*  Find  out  If  the  projection  of  the  * 

*  Sat  1  vector  into  Sat  2*s  * 

*  perifocal  plane  is  within  Sat  2*8  * 

*  orbit.  ♦ 
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ANGLE  «  DATAN2 (XI (2) , XI ( 1) ) 

DIST  *  A2*  ( 1D0-E2*E2)  /  ( 1DO+E2*DCOS  (ANGLD ) 

IF  (DIST  .GT.  PR)  THEN 
INSIDE  -  .TRUE. 

ELSE 

INSIDE  «  .FALSE. 

ENDIF 

IF  (DEBUG)  WRITE (*, 1040)  RADDEG*ANGLE,  DIST,  INSIDE 


♦ - * 

*  If  the  projection  of  Sat  l’s  * 

*  vector  into  Sat  2's  perifocal  * 

*  plane  is  within  Sat  2's  orbit.  * 

*  then  use  the  first  set  of  quadrant  * 

*  checks.  If  outside  of  Sat  2's  * 

*  orbit,  then  use  the  second  set  of  * 

*  quadrant  checks.  * 

*  - * 


NU11  -  0D0 
NU12  -  0D0 
NU21  -  0D0 
NU22  *  0D0 
IF  (INSIDE)  THEN 

IF  (XI (2) .LT.ODO)  THEN 
IF  (PDTH.GT.BOUND1)  THEN 
IF  (PDTH.GT.BOUND2)  THEN 
IF  (N.EG.O)  THEN 
NU1 1  -  0D0 
NU12  -  TWOPI 
ELSE  IF  (N.EQ.2)  THEN 
NU11  -  RANK (2)  -  TWOPI 
NU12  -  RANK(l) 

ELSE  IF  (N.EQ.4)  THEN 
NUU  -  RANK (4)  -  TWOPI 
NU12  ■  RANK(l) 

NU21  -  RANK (2) 

NU22  -  RANK (3) 

ELSE 

WRITER.  1050)  N 
ENDIF 
ELSE 

IF  (N.EQ.2)  THEN 
NU1 1  »  -RANK (2) 

NU12  -  RANK ( 1 ) 

ELSE  IF  (N.EQ.4)  THEN 
NUl 1  -  -  RANK (4) 

NU12  -  RANK(l) 
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NU21  -  RANK (2) 

NU22  -  RANK (3) 

ELSE 

WRITE(* , 1060)  N 
ENDIF 
ENDIF 
ELSE 

IF  (PDTH.GT.B0UND2)  THEN 
IF  (N.EQ.2)  THEN 
NU11  -  RANK (2) 

NU12  -  TWOPI  -  RANK(l) 
ELSE  IF  (N.EQ.4)  THEN 
NU11  -  RANK (2) 

NU12  -  RANK (3) 

NU21  -  RANK (4) 

NU22  -  TWOPI  -  RANK(l) 
ELSE 

WRITEC* , 1070)  N 
ENDIF 
ELSE 

IF  (N.EQ.2)  THEN 

NU11  *  TWOPI  -  RANK (2) 
NU12  -  TWOPI  -  RANK(l) 
ELSE  IF  (N.EQ.4)  THEN 
NU11  *  RANK (2) 

NU12  *  RANK (3) 

NU21  -  TWOPI  -  RANK (4) 
NIJ22  •  TWOPI  -  RANK(l) 
ELSE  IF  (N.NE.O)  THEN 
WRITER.  1080)  N 
ENDIF 
ENDIF 
ENDIF 
ELSE 

IF  (?DTH.GT. BOUND 1)  THEN 
IF  (PDTH . GT . B0UND2)  THEN 
IF  (N.EQ.O)  THEN 
NUil  «  ODO 
NU12  -  TWOPI 
ELSE  IF  (N.EQ.2)  THEN 
NU11  -  -  RANK ( 1 ) 

NU12  *  TWOPI  -  RANK (2) 
ELSE  IF  (N.EQ.4)  THEN 
NU11  «  -  RANK(l) 

NU12  ■  TWOPI  -  RANK (4) 
NU21  -  TWOPI  -  RANK (3) 
NU22  -  TWOPI  -  RANK (2) 
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ELSE 

WRITEC*, 1090)  N 
ENDIF 
ELSE 

IF  CN.EQ.2)  THEN 
NU11  -  -  RANK(l) 

NU12  -  RANK (2) 

ELSE  IF  (N.EQ.4)  THEN 
NU11  -  -  RANK(l) 

NU12  -  RANK (4) 

NU21  -  TWOPI  -  RANK (3) 
NU22  *  TWOPI  -  RANK (2) 
ELSE 

WRITEC*. 1100)  N 
ENDIF 
ENDIF 
ELSE 

IF  (PDTH . GT . B0UND2)  THEN 
IF  CN.EQ.2)  THEN 
NU11  -  RANK(l) 

NU12  *  TWOPI  -  RANKC2) 
ELSE  IF  (N.EQ.4)  THEN 
NUU  -  RANK(l) 

NU12  -  TWOPI  -  RANK (4) 

NU21  ■  TWOPI  -  RANK (3) 

NU22  •  TWOPI  -  RANK (2) 

ELSE 

WRITEC*. 1110)  N 
ENDIF 
ELSE 

IF  CN.EQ.2)  THEN 
NU11  -  RANKC1) 

NU 12  -  RANK (2) 

ELSE  IF  CN.EQ.4)  THEN 
NU 1 1  »  RANKC1) 

NU12  -  RANK (4) 

NU21  -  TWOPI  -  RANK (3) 
NU22  -  TWOPI  •  RANK (2) 
ELSE  IF  CN.NE.O)  THEN 
WRITEC*. 1120)  N 
ENDIF 
ENDIF 
ENDIF 
ENDIF 
ELSE 

IF  (XI (2) .LT.ODO)  THEN 
IF  (PDTH. GT. BOUND 1)  THEN 
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IF  (PDTH.GT.B0UND2)  THEN 
IF  (N.EQ.O)  THEN 
NU11  «  ODO 
NU12  -  TWOPI 
ELSE  IF  (N.EQ.2)  THEN 
NU11  «  HANK (2)  -  TWOPI 
NU12  -  RANKC1) 

ELSE  IF  (N.EQ.4)  THEN 
NU11  -  RANK (4)  -  TWOPI 
NU12  -  RANK(l) 

NU21  -  RANK (2) 

NU22  -  RANK (3) 

ELSE 

WRITER.  1130)  N 
ENDIF 
ELSE 

IF  (N.EQ.2)  THEN 
NU11  -  -  RANK (2) 

NU12  •  RANK(l) 

ELSE  IF  (N.EQ.4)  THEN 
NUU  -  -  RANK (4) 

NU12  -  RANK(l) 

NU21  -  RANK (2) 

NU22  -  RANK (3) 

ELSE 

WRITER.  1140)  N 
ENDIF 
ENDIF 
ELSE 

IF  (PDTH.GT.B0UND2)  THEN 
IF  (N.EQ.2)  THEN 
NU11  -  RANK (2) 

NU12  «  TWOPI  -  RANK(l) 
ELSE  IF  (N.EQ.4)  THEN 
NUU  -  RANK (2) 

NU12  -  RANK (3) 

NU21  -  RANK (4) 

NU22  »  TWOPI  -  RANK(t) 
ELSE 

WRITE(* . 1 150)  N 
ENDIF 
ELSE 

IF  (N.EQ.2)  THEN 

NU11  «  TWOPI  -  RANK (2) 
NU12  -  TWOPI  -  RANK(l) 
ELSE  IF  (N.EQ.4)  THEN 
NUil  -  RANK (2) 
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NU12  -  HANK (3) 

NU21  -  TNOPI  -  RANK (4) 
NU22  «  TNOPI  -  RANK(i) 
ELSE  IF  CN.NE.O)  THEN 
WRITE (*.1160)  N 
ENDIF 
ENDIF 
ENDIF 
ELSE 

IF  (PDTH.GT. BOUND 1)  THEN 
IF  (PDTH.GT. B0UND2)  THEN 
IF  (N.EQ.O)  THEN 
NUU  -  ODC 
NU12  *  TNOPI 
ELSE  IF  (N.E0.2)  THEN 
NU11  -  -  RANKC1) 

NU12  -  TNOFI  -  RANK (2) 
ELSE  IF  (N.E3.4)  THEN 
NU11  «  -  RANK(l) 

NU12  -  TNOPI  -  RANKC4) 
NU21  *  TNOPI  -  RANK (3) 
NU22  «  TNOPI  -  RANK (2) 
ELSE 

WRITER.  1170)  N 
ENDIF 
ELSE 

IF  CN.EQ.2)  THEN 
NUU  «  -  RANKU) 

NU12  -  RANK (2) 

ELSE  IF  (N.EQ.4)  THEN 
NUU  -  -  RANKU) 

NU 12  -  RANK <4) 

NU21  *  TNOPI  -  RANK (3) 
NU22  -  TNOPI  -  RANK (2) 
ELSE 

WRITEC* . 1 180)  N 
ENDIF 
ENDIF 
ELSE 

IF  (PDTH.GT. B0UND2)  THEN 
IF  (N.FQ.2)  THEN 
NUU  -  RANKU) 

NU12  «  TNOPI  -  RANK (2) 
ELSE  IF  (N.E0.4)  THEN 
NU1 1  -  RANKU) 

NU12  -  TNOPI  -  RANK (4) 
NU2!  -  TNOPI  -  RANK (3) 
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NU22  *  TWOPI  -  RANK (2) 

ELSE 

WRITE (♦. 1190)  N 
ENDIF 
ELSE 

IF  (N.EQ.2)  THEN 
NU11  -  RANK(l) 

NU12  *  RANK (2) 

ELSE  IF  (N.EQ.4)  THEN 
NU1 1  -  RANK ( 1 ) 

NU12  -  RANKC4) 

NU21  •  TWOPI  -  RANK (3) 

NU22  -  TNOPI  -  RANK (2) 

ELSE  IF  (N.NE.O)  THEN 
WRITER.  1200)  N 
ENDIF 
ENDIF 
ENDIF 
ENDIF 
ENDIF 
RETURN 

900  FORMAT  <//’  DSQR  *’ .D17. 10.4X. 'PR  ,D17. 10. 4X. ’E2*’ . 

1  D17.10/*  SAT2P-’.D1710.4X.  *X(1>-\D17.I0.4X// 

2  1  A  •’ .D17. 10.4X. *B  *’ .D17. 10.4X. *C  •’ .017.10/ 

3  1  K  *\D17. 10//> 

1000  FORMAT  <//’  Norder  -  Ml.//1  PI** ,D17. 10.4X.  ’P2»* . 

1  D17. 10/'  P3-\D17.10.4X/P4-\D17.10/ 

2  *  P5*’ ,D17. 10//) 

1010  FORMAT  <//’  PC  * . 1 1 . *)»' ,D17. 10.4X. ’ANSWER! ' . II . *)•' . 

1  * ( ' .F16. 10. * . ' .Fl€. 10. *> ’> 

1020  FORMAT  (//■  I-M1//C  Solutlon< ’  .11 .  ’>•*  .F16. 10.  •  deg*/)) 
1030  FORMAT  (//’  N»\II//<’  RankC 1 . II .  *)-  .r  16. 10.  •  dtg1/)) 

1040  FORMAT  <//*  Sc.t  1  Projection  occur*  at  a  Sat  2  NU  of  ’, 

1  F16.10.’  dtg  where  Sat  2  radiu*  It  ’.F16.10.*  ka’// 

2  ’  INSIDE  •  Ml//) 

1050  FORMAT  <//’  ERROR!  N»  * .  1 1 . '  XI  (2)  <  0  PDTH  >  BOUNDl  #. 

1  ’PDTH  >  B0UND2  INSIDE*//) 

1060  FORMAT  (//*  ERROR!  N»\II.’  XI  (2)  <  0  PDTH  >  BOUNDl  \ 

1  *PDTH  <  B0UND2  INSIDE*//) 

1070  FORMAT  <//’  ERROR!  N-’.I!.’  XI  (2)  <  0  PDTH  <  BOUNDl  *. 

1  ’PDTH  >  B0UND2  INSIDE*//) 

1080  FORMAT  (//*  ERROR!  N*‘.I1.’  XI  <2)  <  0  PDTH  <  BOUND!  *. 

1  ’PDTH  <  B0UND2  INSIDE’//) 

1090  FORMAT  <//  ERROR!  N*\Il.‘  X:  C2)  >  0  PDTH  >  BOUND!  \ 

1  ’PDTH  >  B0UND2  INSIDE’//) 

1100  FORMAT  <//*  ERROR!  N-’.Il.’  XI U.  >  0  PDTH  >  BOUNDl  *. 

1  ‘PDTH  <  BOUND2  INSIDE*//) 


o  o  o  o  o 


1110  FORMAT  (//’  ERROR  I  N-Ml,’  XI  (2)  >  0 
1  TDTH  >  BOUND2  INSIDE’//) 

1120  FORMAT  (//’  ERROR!  N-’.Il/  XI  (2)  >  0 
1  ’PDTH  <  BQUND2  INSIDE’//) 

1130  FORMAT  (//’  ERROR!  N-Ml.’  XI  (2)  <  0 
1  ’PDTH  >  B0UND2  OUTSIDE’//) 

1140  FORMAT  <//’  ERROR!  N-Ml.’  XI (2)  <  0 
1  TDTH  <  BOUND2  OUTSIDE’//) 

1150  FORMAT  <//’  ERROR!  N-’.I1.'  XI  (2)  <  0 
1  TDTH  >  B0UND2  OUTSIDE’//) 

1160  FORMAT  (//’  ERROR!  N-Ml.’  XI  (2)  <  0 
1  TDTH  <  B0UND2  OUTSIDE’//) 

1170  FORMAT  (//’  ERROR!  N«’.I1.’  XI (2)  >  0 
1  TDTH  >  B0UND2  OUTSIDE’//) 

1180  FORMAT  <//’  ERROR!  N-Ml.’  XI  (2)  >  0 
1  TDTH  <  B0UND2  OUTSIDE’//) 

1190  FORMAT  (//’  ERROR!  N-Ml.’  XI  (2)  >  0 
1  TDTH  >  B0UND2  OUTSIDE*//) 

1200  FORMAT  (//’  ERROR!  N-Ml.’  XI  (2)  >  0 
1  ’PDTH  <  B0UND2  OUTSIDE’//) 

END 


PDTH  <  BOUND 1  ’ 
PDTH  <  BOUND 1  ’ 
PDTH  >  BOUND 1  ’ 
PDTH  >  BOUND 1  ’ 
PDTH  <  BOUND 1  ’ 
PDTH  <  BOUND1  ’ 
PDTH  >  BOUND 1  ’ 
PDTH  >  BOUND 1  ’ 
PDTH  <  BOUND 1  ’ 
PDTH  <  BOUND 1  ’ 


Function  MO 


REAL*8  FUNCTION  MCNU) 

C  ♦ . . 

10GICAL*4  DEBUG 

REAL'S  NU.  SHNU.  CHNU.  X.  Y.  Z 

REAL'S  BOUND 1 .  B0UND2.  DTH.  DSQR.  PDTH.  PR.  IU11.  NU12. 
1  NU21.  NU22 

REAL'8  Ml.  Al.  El.  INCi .  NU1 .  ARGPA1 .  LONAN1 .  XI (4). 

1  SINC1.  CINC1.  SNU1 .  CNUl .  SARGP1 .  CARGP1 . 

2  SLNANi.  CLNAN1.  BETA! 

REALMS  M2.  A2.  E2.  INC2.  NU2.  ARGPA2.  IONAH2.  X2(4) . 

1  SINC2.  CINC2.  SNU2.  CNU2.  SARGP2.  CARCP2. 

2  SLNAN2.  CLNAN2.  BETA2 

REAL'8  Rll.  R12.  Ri3.  R21 .  R22.  R23.  R31 .  R32.  R33. 

1  Sll.  S12.  S13.  S21.  S22.  S23.  S31 .  S32.  S33 
REAL'S  HAIFPI . PI . TNOPI . DEGRAD . RADDEG 
REAL'S  ER.  DU.  TU 
COMMON  /ADMIN/  DEBUG 


COMMON  /ASTRO/ 
COMMON  /BOUND/ 


ER.  DU.  TU 

BOUND 1 .  B0UND2.  DTH.  DSQR.  PDTH.  PR.  NUll. 


1  NU 12.  NU21 .  NU22 

COMMON  /CONST/  HALF? I .  PI.  TWOPI .  DEGRAD.  RADDEG 


1U 


o  o  o  o  o 


COMMON  /SAT1/  Ml.  A1 ,  El.  INC1 .  NU1 .  ARGPA1 ,  LONAN1 .  XI, 

1  SINC1 ,  CINC1 ,  SNU1 ,  CNU1 ,  SARGP1 .  CARGP1 . 

2  SLNAN1 .  CLNAN1 ,  BETA1 

COMMON  /SAT2/  M2.  A2.  E2,  INC2,  NU2,  ARGPA2 ,  L0NAN2 ,  X2. 

1  SINC2,  CINC2,  SNU2,  CNU2,  SARGP2 ,  CARGP2 . 

2  SLNAN2 .  CLNAN2 ,  BETA2 

COMMON  /TRANS/  Rll,  R12,  R13,  R21 .  R22.  R23.  R31 ,  R32. 

1  R33,  SI  1 .  S12.  S13,  S21 ,  S22.  S23.  S31 .  S32.  S33 

C  * . * 

IF  (NU.EQ.ODO)  THEN 
M  -  0D0 
RETURN 

ELSE  IF  (DABS(NU) .EQ.PI)  THEN 
M  -  DSIGN(PI.NU) 

RETURN 

ELSE  IF  (NU.EQ.TNOPI)  THEN 
M  «  TWOPI 
RETURN 
ENDIF 

SHNU  -  DSINC5D-1+NU) 

CHNU  -  DCOS(5D-l*NU) 

IF  (CHNU.EQ.ODO)  THEN 
M  -  DSIGN(PI.NU) 

IF  (DEBUG)  WRITE(* , 1000)  RADDEG*NU ,  RADDEG*M 
RETURN 
ENDIF 

M  *  2D0*(  DATANC  BETA2*SHNU/CHNU  )  -  E2*BETA2*SHNU*CHNU/ 

1  (CHNU*CHNU  ♦  BETA2*BETA2  *  SHNU*SHNU  )) 

IF  (M.LT.ODO)  M  -  TNOPI  ♦  M 

IF  (NU.LT.ODO)  M  ■  M  -  TNOPI 

IF  (DEBUG)  NRITE(* . 1000)  RADDEG*NU .  RADDEG*M 

RETURN 

1000  FORMAT  ('  NU- 1  .F16 . 10 . 1  deg\4X.'  M-\F16.10.'  deg*) 

END 


Function  PDCAO 


REAL*8  FUNCTION  PDCA(NU) 

C  * . * 

L0GICAL*4  DEBUG 
INTEGER*4  I 

REAL*8  M.  SHNU  1 ,  CHNUl.  DMPDN’J 

REAL*8  NU.  DX.  DY.  X.  Y.  Z 

REAL*8  B0UND1 .  B0UND2.  DTH,  DSQR.  PDTH.  PR.  NU11.  NU12. 
1  NU21.  NU22 


113 


o  o  o  o  oooo  ooo 


REAL*8  Ml.  Al.  El.  INC1 ,  NU1 .  ARGPA1 .  L0NAN1,  XI (4), 

1  SINC1 .  CINC1 ,  SNU1 ,  CNU1 ,  SARGP1 ,  CARGP1 . 

2  SLNAN1 ,  CLNAN1,  BETA1 

REAL*8  M2,  A2,  E2,  INC2,  NU2,  ARGPA2,  L0NAN2 ,  X2(4) , 

1  SINC2 .  CINC2 ,  SNU2,  CNU2,  SARGP2 ,  CARGP2 , 

2  SLNAN2,  CLNAN2,  BETA2 

REAL*8  Rll,  R12,  R13,  R21 .  R22,  R23,  R31 ,  R32.  R33, 

1  Sll.  S12,  S13.  S21,  S22,  S23,  S31 ,  S32,  S33 
REAL*8  HALFPI . PI , TNOPI . DEGRAD , RADDEG 
REAL*8  ER,  DU,  TU 
COMMON  /ADMIN/  DEBUG 
COMMON  /ASTRO/  ER,  DU,  TU 

COMMON  /BOUND/  BOUND 1 .  B0UND2 ,  DTH,  DSQR,  PDTH,  PR.  NU11, 
1  NU12  NU21  NU22 

COMMON* /CONST/  HALFPI,  PI.  TWOPI ,  DEGRAD,  RADDEG 
COMMON  /SAT 1 /  Ml,  Al .  El.  INC1 ,  NU1 .  ARGPA1 ,  LONAN1 .  XI. 

1  SINC1.  CINC1 ,  SNU1 ,  CNU1 ,  SARGP1 ,  CARGP1 , 

2  SLNAN1 ,  CLNAN 1 .  BETA1 

COMMON  /SAT2/  M2.  A2.  E2,  INC2.  NU2,  ARGPA2,  L0NAN2 ,  X2. 

1  SINC2,  CINC2,  SNU2,  CNU2,  SARGP2 .  CARGP2, 

2  SLNAN2 ,  CLNAN2 ,  BETA2 

COMMON  /TRANS/  Rll,  R12,  R13.  R21 ,  R22,  R23,  R31 ,  R32, 

1  R33,  Sll,  S12,  S13.  S21,  S22.  S23,  S31 .  S32,  S33 


C  * - - * 

NU1  -  NU 
SNU1  -  DSIN(NUl) 

CNU1  -  DCOS(NUl) 

* - — - * 

*  Compute  the  Radius  of  Sat  1  * 

* — - - - - * 

XI C4)  -  Al  *  (IDO  -  E1*E1)  /  (IDO  ♦  E1*CNU1) 

*  - - - * 

*  Compute  the  position  of  Sat  1  in  * 

*  Perifocal  Coordinate  Frame  * 

*  - — - — - - — — * 

XI (1)  -  XI (4)*CNU1 
XI (2)  -  XI (4) *SNU1 
XI (3)  -  0D0 


IF  (DEBUG)  WRITER,  1000)  RADDEG*NU1 .  (XI  (I)  .1-1 .4) 


* - — „ — - — - — - * 

*  Compute  the  position  of  Sat  1  in  * 

*  the  Geocentric  Equatorial  Frame.  * 

*  - - - * 

X  -  R1 1*X1 (1)  ♦  R12*X1 (2) 

Y  -  R21*X1 (1)  ♦  R22*X1 (2) 

Z  -  R31*X1 (1)  ♦  R32*X1 (2) 


IF  (DEBUG)  ¥RITE(* . 1010)  X . Y ,Z ,DSQRT(X*X*Y*Y*Z*Z) 
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C  * - * 

C  *  Transform  Sat  l's  position  from  * 
C  *  the  Geocentric  Equatorial  Frame  * 
C  *  to  Sat  2*s  perifocal  frame  * 

C  * - " - * 

XI (1)  «  S11*X  ♦  S12*Y  ♦  S13*Z 
XI (2)  -  S21*X  ♦  S22*Y  ♦  S23*Z 
XI (3)  *  S31*X  ♦  S32*Y  ♦  S33*Z 


PR  -  DSQRT(X1(1)*X1  (1)4X1  (2)*X1  (2)) 

IF  (DEBUG)  WRITEC*, 1020)  CXI  Cl) .1-1 .3) . 

1  DSQRTCX1 (1)*X1 (1)*X1 (2)*X1 (2)*X1 (3)*X1 (3)).  PR 


* - - — - - * 

*  Compute  the  distance  between  XI  * 

*  projected  onto  Sat  2’s  perifocal  * 

*  plane  and  Sat  2’s  perigee  (BOUND1)  * 

*  and  apogee  (B0UND2)  * 

*  - - — - - - — - — — — * 


DX  -  XI (1)  -  A2* (1D0-E2) 

DY  -  XI (2) 

BOUND1  -  DSQRT (DX*DX  ♦  DY*DY) 
DX  -  XI (1)  ♦  A2*(1D0*E2) 
B0UND2  -  DSQRT (DX*DX  ♦  DY*DY) 


*  Find  out  if  there  are  any  points 

*  on  Sat  2‘s  orbit  that  are  exactly 

*  DTH  away  from  the  endpoint  of 

*  vector  XI.  If  there  are.  then 

*  there  are  either  two  or  four 

*  points. 


DSQR  •  DTH*DTH  -  XI (3) *X 1(3) 

PDTH  •  0D0 

IF  (DSQR  .GT.  0D0)  THEN 
PDTH  «  DSQRT (DSQR) 

IF  (DEBUG)  MRITE(* , 1030)  PDTH.  BOUND 1 .  BOUND2 
CALL  ORBBNDO 
IF  (DEBUG)  CALL  PRINTRO 
ELSE 

IF  (DEBUG)  WRITER,  1030)  PDTH.  BOUND  1 .  BOUND2 


PDCA  -  0D0 
RETURN 
ENDIF 

* - — - - ——————— — - — -* 

♦  Compute  Probability  Density  of  * 

♦  Close  Approach  * 

♦  - - - — - - * 


PDCA  -  (M(NU12) -M(NU1 1)  ♦  M(NU22)-M(NU21))  /(TWOPI*TNOPI) 
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IF  (DEBUG)  WRITE(* .  1040)  PDCA 
SHNU1  -  DSIN (5D-1*NU1) 
CHNU1  -  DC0S(5D-1*NU1) 


DMPDNU  -  CHNU 1 *CHNU 1  ♦  BETA 1* BETA 1*SHNU1*SHNU1 
DMPDNU  «  ( 1 D0-E1 ) *BETA 1 / (DMPDNU*DMPDNU) 

PDCA  «  PDCA*DMPDNU 

IF  (DEBUG)  WRITE (* , 1050)  PDCA,  DMPDNU 
RETURN 

1000  FORMAT  (//’  Input  NU«\F16.10// 

1  ’  Sat  1  Position  in  Sat  l”s  perifocal  frame:*// 

2  *  X  Axis:  \F17.10,’  km’/ 

3  ’  Y  Axis:  ’.F17.10,  *  km*/ 

4  ’  Z  Axis:  ’.F17.10.  *  km’/ 

5  ’  Sat  1  Radius:  ’.F17.10.’  km’//) 


1010  FORMAT  (’  Sat  1  Position  in  the  inertial  reference  frame 


1 

2 

3 

4 


//’ 


Y  Axis 
Z  Axis 
Sat  1  Radius 


X  Axis:  * ,F17 .  10 .  *  km’/ 
’ . F17 . 10 . ’  km’/ 

’ , F17. 10 . ’  km’/ 

’ , F17 . 10 . *  km’//) 


1020  FORMAT  (’  Sat  1  Position  in  Sat  2’’s  Derifocal  frame:’// 

1  *  X  Axis:  ’ . F17 . 10 . ’  km’/ 

2  ’  Y  Axis:  * . F17 . 10 , ’  km’/ 

3  ’  Z  Axis:  ’.F17.10.’  km’/ 

4  ’  Sat  1  Radius:  ’ .F17. 10. ’  km’/ 

5  ’  Projection:  ’.F17.10,*  km’/ 

6  ’  Radius’//) 

1030  FORMAT  (//’  PDTH-’ .F17. 10/’  BOUND1-* ,F17. 10. 


1  4X.*BOUND2-’.F17.10/) 

1040  FORMAT  (’  Unsealed  PDCA-’ . FI 7. 10/) 

1050  FORMAT  (’  PDCA-’ .F16. 10. 4X. ’DMPDNU-’ .F16. 10//) 
END 


C 
C 
C 
C 
C 

SUBROUTINE  PRINTRO 


Subroutine  PRINTRO 


C  ♦ . * 

LOG I CAL* 4  DEBUG 
INTEGER*4  I.  J.  K 

REAL*8  DX.  DY.  DZ.  DR1 .  DR2.  DR3.  DR4 
REAL*8  BOUND1 .  BOUND2.  DTH.  DSQR.  PDTH,  PR.  NUll,  NU12. 
1  NU21 .  NU22 


REAL*6  Ml.  Al.  El.  INCl .  NU1 .  ARGPA1 .  LQNANl .  XI (4) 
REAL*8  M2.  A2.  E2.  INC2.  NU2,  ARGPA2.  L0NAN2.  X2(4) 
REAL*8  SINC1 .  CINCl .  SNU1.  CNUl.  SARGPl .  CARGP1 . 

1  SLNAN1 .  CLNAN1.  BETAl 
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REAL*8  SINC2,  CINC2,  SNU2,  CNU2,  SARGP2 ,  CARGP2 , 

1  SLNAN2 ,  CLNAN2 ,  BETA2 
REAL+8  HALFPI . PI , TWOPI , DEJRAD , RADDEG 
REAL+8  ER,  DU,  TU 
COMMON  /ADMIN/  DEBUG 
COMMON  /ASTRO/  ER,  DU,  TU 

COMMON  /BOUND/  BOUND 1 .  B0UND2 ,  DTH,  DSQR,  PDTH,  PR,  NU1 1 , 
1  NU12,  NU21,  NU22 

COMMON  /CONST/  HALFPI.  PI,  TNOPI ,  DEGRAD,  RADDEG 
COMMON  /SAT1/  Ml.  A1 .  El,  INC1 ,  NU1 .  ARGPA1 .  LONAN1 ,  XI, 

1  SINC1,  CINC1.  SNU1 .  CNU1 ,  SARGP1 ,  CARGP1 , 

2  SLNAN1.  CLNAN 1 ,  BETA1 

COMMON  /SAT2/  M2.  A2,  E2,  INC2,  NU2,  ARGPA2 ,  L0NAN2,  X2, 

1  SINC2,  CINC2,  SNU2.  CNU2,  SARGP2 ,  CARGP2, 

2  SLNAN2 .  CLNAN2 .  BETA2 


X2(4)  -  A2*(1DO-E2*E2)/(1DO+E2*DCOS(NU1 i)> 

X2(l)  «  X2(4) *DCOS(NUl 1) 

X2C2)  -  X2(4)*DSIN (NU1 1) 

X2C3)  •  0D0 

DX  -  X2C1)  -  XI (1) 

DY  -  X2C2)  -  XI (2) 

DZ  •  X2(3)  -  XI (3) 

DR1  «  DSQRT(DX*DX  ♦  DY*DY  ♦  DZ*DZ) 

X2(4)  -  A2*(1DO-E2*E2)/(1DO+E2*DCOS(NU12)> 

X2(l)  -  X2C4)*DCOS(NU12) 

X2(2)  -  X2(4)*DSIN(NU12) 

X2(3)  -  0D0 

DX  •  X2(l>  -  XI <1) 

DY  »  X2(2)  -  XI (2) 

DZ  •  X2(3)  -  XI (3) 

DR2  -  DSQRT <DX*DX  ♦  DY*DY  ♦  DZ*DZ) 

WRITEC* . 1000)  RADDEG*NU1 1 ,  DR1 .  RADDEG*NU12.  DR2 
X2(4)  -  A2*(1D0-E2*E2)/(1D0*E2*DCQS(NU21)> 

X2U)  *  X2(4)*DC0S(NU21) 

X2(2)  •  X2(4)*DSIN(NU2l) 

X2C3)  -  ODO 

DX  «  X2<1)  -  XI (1) 

DY  -  X2(2)  -  XI (2) 

DZ  -  X2<3)  -  XI (3) 

DR3  -  DSQRT (DX*DX  ♦  DY*DY  ♦  DZ*DZ> 

X2<4)  -  A2*UD0-E2*E2>/C1D0*E2*DCQS(NU22>> 

X2<1)  -  X2<4>*DCGS(NU22) 

X2C2)  «  X2(4) *DSIN  CNU22) 

X2C3)  -  ODO 

DX  -  X2U)  -  XI  (1) 

DY  -  X2(2)  -  XI <2) 
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DZ  *  X2(3)  -  XI (3) 

DR4  »  DSQRT <DX*DX  ♦  DY*DY  ♦  DZ*DZ) 

WRITE<*,1100)  RADDEG*NU21 ,  DR3,  RADDEG*NU22 ,  DR4 
RETURN 

1000  FORMAT  (//’  NU11=’  .F16. 10.  ’  deg’.4X. 

1  ’Range  to  NU11  *  ’.F16.10,*  km’/ 

2  ’  NU12*’ ,F16. 10, ’  deg* ,4X, 

3  ’Range  to  NU12  «  ’.F16.10,’  km’) 

1100  FORMAT  (’  NU21*’ ,F16. 10, ’  deg’,4X. 

1  ’Range  to  NU21  «  \F16.10.’  km’/ 

2  ’  NU22*’ ,F16 . 10 , , ’  deg’4X. 

3  ’Range  to  NU22  «  ’.F16.10.’  km’//) 

END 


*  Function  CCUBRT  * 

*  * 

* - - - — * 

COMPLEX* 16  FUNCTION  CCUBRT (X) 

COMPLEX* 16  X 

C  * . * 

REAL*8  A.  ANGLE.  B.  MAG 
REAL*8  HALFPI . PI . TWOPI . DEGRAD . RADDEG 
COMMON  /CONST/  HALFPI.  PI ,  TWOPI.  DEGRAD.  RADDEG 
C  * . . . * 


A  «  DREAL(X) 

B  «  DIMAG (X) 

IF  (  B  ,EQ.  0D0)  THEN 
MAG  -  DABS (A) 

MAG  «  DSIGN<  MAG**(1D0/3D0) .  A) 

CCUBRT  «  DCMPLXC  MAG.  0D0) 

ELSE 

MAG  «  CDABS(X) 

ANGLE  *  DATAN2(  B.  A) 

IF  (A  .LT.  0D0)  THEN 
ANGLE  »  CDSIGN (TWOPI .ANGLE) *ANGLE) /3D0 
ELSE 

ANGLE  «  ANGLE/3D0 
ENDIF 

MAG  «  MAG**(1D0/3D0) 

CCUBRT  «  DCMPLX CMAG*DCOS (ANGLE)  .MAG*DSIN (ANGLE)) 
ENDIF 
RETURN 
END 


C  * . . . - . • 

c  •  • 

C  *  Subroutine  QUADRADIC  * 
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c  *  * 

C  * - * 

SUBROUTINE  QUADRADIC(P.X) 

REAL+8  P<2> 

COMPLEX* 1 6  X(2) 

C  * - * 

REAL*8  A 
COMPLEX* 16  B 

C  * - * 

A  -  -PCD 


B  -  CDSQRT (DCMPLX (A*A  -  4D0*P(2))) 

X(l)  -  (A  -  B)/2D0 

XC2)  -  (A  ♦  B)/2D0 

RETURN 

END 


* - - - — * 

*  * 

*  Subroutine  CUBIC  * 

♦  * 

* - — - - - — * 

SUBROUTINE  CUBIC (P.X) 

REAL*8  PC3) 

COMPLEX* 16  X(3) 

C  * . * 

REAL*8  A.  B.  C 

COMPLEX* 16  CCUBRT .  CX.  D.  E.  SQRT3XJ 
C  * . * 


SQRT3XJ  *  (ODO.  1 .73205080756887729DO) 

A  -  (3D0*P(2)~P(1)*P(1))/3D0 

B  «  (2D0*PU)*PU)*P(i)  -  9D0*P(1)*P(2)  ♦  27D0*P(3))/27D0 
C  -  B*B/4DO  ♦  A*A*A/27D0 
IF  (C.LT.ODO)  THEN 
CX  *  (ODO . 1 DO) *DSQRT (DABS (C) ) 

ELSE 

CX  «  (IDO. ODO) *DSQRT (C) 

ENDIF 

D  «  CCUBRT (-B/2D0  ♦  CX) 

E  «  CCUBRT (-B/2D0  -  CX) 

X(l)  «  D  ♦  E  -  P(1)/3D0 

X(2)  -  -(D  ♦  E)/2D0  ♦  (D  -  D/2D0*SQRT3XJ  -  P(1)/3D0 
X(3)  -  -(D  ♦  E)/2D0  -  (D  -  E)/2D0*S0RT3XJ  -  P(1)/3D0 
RETURN 
END 


C  ♦ . - . * 

C  *  * 

C  *  Subroutine  OUART1C  * 

C  *  * 
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SUBROUTINE  QUARTIC (P . X) 

REAL*8  P(4) 

COMPLEX* 16  X(4) 

C  * - * 

INTEGER  I.  IS 
REAL*8  CC3) .  Y 
COMPLEX* 16  D.  E.  R.  SC3) 

C  * - * 

C(l)  -  -P(2) 

C(2)  -  P(1)*P(3)-4D0*PC4> 

C(3)  *  -P(l) *P(1) *P(4)  ♦  4D0*P(2)*P(4)  -  P(3)*P(3) 
CALL  CUBIC (C.S) 

IS  -  0 

DO  100  1-1.3 

IF  (DIMAG (S(I)).EQ.0D0)  THEN 
IF  CIS. EQ. ODO)  THEN 
IS  »  I 

Y  -  DREAL(SCI)) 

ELSE 

IF  (DREALCSCD).GT.Y)  THEN 
IS  »  I 

Y  -  DREAL (SCI)) 

ENDIF 

ENDIF 

ENDIF 

100  CONTINUE 

IF  (IS.EO.O)  THEN 
Y  ■  ODO 
NRITEC* . 1000) 

DO  200  >1.3 

200  WRITEC* . 1010)  DREAL (S(I)) .  DIMAG (SCI)) 

ENDIF 

R  •  CDSQRT (DCMPLX (P ( 1 ) *P ( 1 ) /4D0  -  P(2)  ♦  Y)) 

IF  (R.EQ. (ODO. ODO))  THEN 
D  -  2D0*CDSQRT (DCMPLX ( Y* Y  -  4D0*P(4))) 

E  -  -D 

D  -  CDSQRT ( . 75D0*P ( 1 ) *P ( 1 )  -  2D0*P(2)  ♦  D) 

E  -  CDSQRT ( .75D0*PU)  *P(l)  -  2D0*P(2)  ♦  E) 

ELSE 

D  •  (4D0*PU)*PC2)-8D0*P(3)-PU)*P(1)*PU))/4D0/R 
E  »  -D 

D  -  CDSQRT ( .75D0*P(l)*P(l)  -  R*R  -  2D0*P(2>  ♦  D) 

E  -  CDSQRT C.75D0*P(1)*PC1)  -  R*R  -  2D0*PC2)  ♦  E) 
ENDIF 

X(l)  »  -PCD/4D0  ♦  R/2D0  ♦  D/2D0 
X(2)  -  -PCD/4D0  ♦  R/2D0  -  D/2D0 
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XC3)  «  -PCD/4D0  -  R/2D0  ♦  E/2D0 
X(4)  -  -P<1)/4D0  -  R/2D0  -  E/2D0 
RETURN 

1000  FORMAT  (*  Cubic  Error!  All  3  roots  were  comolexl*/) 
1010  FORMAT  (’  Root:  Ml.’  ( *  ,F17. 10.  * .  *  .F17. 10.  * )  •> 

END 
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each  other  at  a  random  time  within  a  specified  time  Interval. 

In  this  paper,  methods  were  developed  to  calculate  probability 
of  close  approach  between  two  satellites.  -To  simplify  the 
analysls.cghe  Investigation  was  restricted  to  satellite  orbits 
and  time  Intervals  where  the  mean  anomaly  of  both  satellites  can 
be  treated  as. Independent,  uniformly  distributed  random 
variables.  SkHoddftloh^U  orbital  parameters,  except  for  mean 
anomaly,  were  assumed  to  be  constant  over  time*.  Thlt- means  that 
all  the  methods  developed  papefrto  calculate  the 

probability  of  close  approach  will  only  be  valid  over  very  long 
time  Intervals  where  the  ratio  of  the  orbital  periods  of  the  two 
satellites  can  be  approximated  as  an  Irrational  number. 

Likewise,  there  can  be  no  perturbations  In  the  orbital 
parameters  of  both  satellites.  , 
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The  first  method  developed  was  $  general  method/for 
calculating  the  probability  of  close  approach  between  two 
satellites  In  elliptical  orbits. t  TheneethodTrequi res  numerical 
integration  and  direct  solution  of  the  roots  of  a  4th  order 
polynomial  during  each  numerical  Integration  step. 

-  -Another  method  was  developed  for  calculating  the  . 
probability  of  close  approach  between  two  satellites  In  circular 
orbits.  This  method  still  requires  numerical  Integration  to 
obtain  a  solution,  but  In  this  case  a  direct  solution  was  found 
for  the  limits  of  Integration.  Furthermore,  the  calculations 
required  during  each  numerical  Integration  step  are  much  simpler 
than  those  required  to  calculate  the  probability  of  close 
approach  with  elliptical  orbits. 

Finally,  a  direct  solution  for  approximate  probability  of 
close  approach  between  two  satellites  In  circular  orbits  was 
developed  for  the  case  where  the  angle  between  the  orbital 
planes  of  both  satellites  Is  not  small  and  the  probability  of 
close  approach  Is  small. 
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s'  Both  the  elliptical  orbit  and  the  circular  orbit  methods  of 
computing  probability  of  close  approach  yielded  results  that 
compare  favorably  with  estimates  of  probability  of  close 
approach  derived  from  statistical  simulations.  ' 


